Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16885) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 16886) @@ -127,8 +127,144 @@ }/*}}}*/ ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/ - _error_("S"); + + /*Intermediaries */ + bool incomplete_adjoint; + IssmDouble Jdet,mu_prime; + IssmDouble eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij; + IssmDouble eps1[3],eps2[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ + element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(element); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + + /*Allocate dbasis*/ + IssmDouble* dbasis = xNew(3*numnodes); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2.*epsilon[1]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; + + for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(dbasis); + xDelete(xyz_list); + return Ke; }/*}}}*/ ElementMatrix* AdjointHorizAnalysis::CreateKMatrixFS(Element* element){/*{{{*/ - _error_("S"); + + /*Intermediaries */ + bool incomplete_adjoint; + IssmDouble Jdet,mu_prime; + IssmDouble eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij; + IssmDouble eps1[3],eps2[3],eps3[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + int numdof = vnumnodes*3 + pnumnodes; + + /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/ + element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(element); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input); + Input* vz_input = element->GetInput(VzEnum);_assert_(vz_input); + + /*Allocate dbasis*/ + IssmDouble* dbasis = xNew(3*vnumnodes); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; + + for(int i=0;ivalues[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; + + Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; + + Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; + Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,XYZEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(dbasis); + xDelete(xyz_list); + return Ke; }/*}}}*/ ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -156,5 +292,5 @@ /*Intermediaries */ int num_responses,i,meshtype,dim; - IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; + IssmDouble Jdet,obs_velocity_mag,velocity_mag; IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; IssmDouble scalex,scaley,scale,S; @@ -355,5 +491,5 @@ /*Intermediaries */ int num_responses,i; - IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; + IssmDouble Jdet,obs_velocity_mag,velocity_mag; IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; IssmDouble scalex,scaley,scale,S; @@ -553,5 +689,5 @@ /*Intermediaries */ int num_responses,i; - IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; + IssmDouble Jdet,obs_velocity_mag,velocity_mag; IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; IssmDouble scalex,scaley,scale,S; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16885) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16886) @@ -172,5 +172,9 @@ virtual void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; + virtual void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; + virtual void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; virtual void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; virtual int VelocityInterpolation()=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16885) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16886) @@ -275,5 +275,5 @@ /*Compute strain rate viscosity and pressure: */ - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -341,5 +341,5 @@ /*Compute strain rate viscosity and pressure: */ - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3d(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -1578,58 +1578,4 @@ return tau_parameter; -} -/*}}}*/ -/*FUNCTION Penta::GetStrainRate3dHO{{{*/ -void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input){ - /*Compute the 3d Blatter/HOStrain Rate (5 components): - * - * epsilon=[exx eyy exy exz eyz] - * - * with exz=1/2 du/dz - * eyz=1/2 dv/dz - * - * the contribution of vz is neglected - */ - - int i; - IssmDouble epsilonvx[5]; - IssmDouble epsilonvy[5]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; -} -/*}}}*/ -/*FUNCTION Penta::GetStrainRate3d{{{*/ -void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ - /*Compute the 3d Strain Rate (6 components): - * - * epsilon=[exx eyy ezz exy exz eyz] - */ - - IssmDouble epsilonvx[6]; - IssmDouble epsilonvy[6]; - IssmDouble epsilonvz[6]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !vz_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); - vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); - - /*Sum all contributions*/ - for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; } /*}}}*/ @@ -3551,5 +3497,5 @@ thickness_input->GetInputValue(&thickness,gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -3574,5 +3520,5 @@ _assert_(gauss->Enum()==GaussPentaEnum); - this->GetStrainRate3d(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); GetPhi(&phi,&epsilon[0],viscosity); @@ -3589,5 +3535,5 @@ IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - this->GetStrainRate3d(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity, &epsilon[0]); @@ -3603,9 +3549,73 @@ IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - this->GetStrainRate3dHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); material->GetViscosity3d(&viscosity, &epsilon[0]); /*Assign output pointer*/ *pviscosity=viscosity; +} +/*}}}*/ +/*FUNCTION Penta::ViscosityHODerivativeEpsSquare{{{*/ +void Penta::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +} +/*}}}*/ +/*FUNCTION Penta::StrainRateHO{{{*/ +void Penta::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){ + /*Compute the 3d Blatter/HOStrain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + int i; + IssmDouble epsilonvx[5]; + IssmDouble epsilonvy[5]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; +} +/*}}}*/ +/*FUNCTION Penta::ViscosityFSDerivativeEpsSquare{{{*/ +void Penta::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +} +/*}}}*/ +/*FUNCTION Penta::StrainRateFS{{{*/ +void Penta::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ + /*Compute the 3d Strain Rate (6 components): + * + * epsilon=[exx eyy ezz exy exz eyz] + */ + + IssmDouble epsilonvx[6]; + IssmDouble epsilonvy[6]; + IssmDouble epsilonvz[6]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vz_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); + vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); + + /*Sum all contributions*/ + for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; } /*}}}*/ @@ -4450,5 +4460,5 @@ GetNodalFunctionsP1(&L[0], gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -4720,5 +4730,5 @@ GetNodalFunctionsP1(&L[0], gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -5361,5 +5371,5 @@ GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; @@ -5432,5 +5442,5 @@ GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; @@ -6845,6 +6855,6 @@ tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); material->GetViscosity3d(&viscosity, &epsilon[0]); material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); @@ -7052,5 +7062,5 @@ GetBprimeSSAFS(&Bprime2[0][0], &xyz_list[0][0], gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity, &epsilon[0]); @@ -7169,5 +7179,5 @@ GetLprimeFSSSA(&LprimeFSSSA[0][0], &xyz_list[0][0], gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -7467,6 +7477,6 @@ if(approximation==SSAHOApproximationEnum){ - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); material->GetViscosity3d(&viscosity, &epsilon[0]); material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); @@ -7475,5 +7485,5 @@ } else if (approximation==SSAFSApproximationEnum){ - this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&newviscosity,&epsilons[0]); } @@ -7710,6 +7720,6 @@ GetBprimeHO(&Bprime[0], &xyz_list[0][0], gauss); - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); material->GetViscosity3d(&viscosity, &epsilon[0]); material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); @@ -7921,5 +7931,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8019,5 +8029,5 @@ GetBprimeFS(Bprime,&xyz_list[0][0],gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8277,5 +8287,5 @@ vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8355,5 +8365,5 @@ NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); @@ -8439,5 +8449,5 @@ vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8520,5 +8530,5 @@ NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); @@ -9042,5 +9052,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -9446,5 +9456,5 @@ GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; @@ -9523,10 +9533,10 @@ GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss); - //this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + //this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]); //eps1[0]=epsilon[0]; eps2[0]=epsilon[3]; eps3[0]=epsilon[4]; //eps1[1]=epsilon[3]; eps2[1]=epsilon[1]; eps3[1]=epsilon[5]; //eps1[2]=epsilon[4]; eps2[2]=epsilon[5]; eps3[2]=epsilon[2]; - this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16885) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16886) @@ -298,5 +298,9 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; + void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); + void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16885) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16886) @@ -173,5 +173,9 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; + void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; + void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16885) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16886) @@ -330,5 +330,9 @@ void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; + void ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; + void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);