Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17545) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17546) @@ -206,4 +206,5 @@ ./shared/Numerics/isnan.cpp\ ./shared/Numerics/cubic.cpp\ + ./shared/Numerics/NewtonSolveDnorm.cpp\ ./shared/Numerics/extrema.cpp\ ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17545) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17546) @@ -3446,4 +3446,5 @@ sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input); } + gauss = element->NewGauss(); for(int i=0;iFindParam(&meshtype,MeshTypeEnum); switch(meshtype){ - case Mesh2DverticalEnum: dim = 2; break; - case Mesh3DEnum: dim = 3; break; - case Mesh3DtetrasEnum: dim = 3; break; + case Mesh2DverticalEnum: dim = 2; tausize = 3; break; + case Mesh3DEnum: dim = 3; tausize = 6; break; + case Mesh3DtetrasEnum: dim = 3; tausize = 6; break; default: _error_("mesh "<Size();i++){ @@ -4265,10 +4274,126 @@ /*Get inputs and parameters*/ - Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + element->GetVerticesCoordinates(&xyz_list); + Input* B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); Input* vz_input; if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} - _error_("NOT implemented yet"); + /*Fetch number of nodes and dof for this finite element*/ + int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG + + /*Initialize vectors*/ + IssmDouble* tbasis = xNew(tnumnodes); + IssmDouble* Ke = xNewZeroInit(tnumnodes*tnumnodes); + IssmDouble* pe_xx = xNewZeroInit(tnumnodes); + IssmDouble* pe_yy = xNewZeroInit(tnumnodes); + IssmDouble* pe_xy = xNewZeroInit(tnumnodes); + IssmDouble* pe_zz = NULL; IssmDouble* pe_xz = NULL; IssmDouble* pe_yz = NULL; + if(dim==3){ + pe_zz = xNewZeroInit(tnumnodes); + pe_xz = xNewZeroInit(tnumnodes); + pe_yz = xNewZeroInit(tnumnodes); + } + + /*Get tau*/ + Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input); + Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input); + Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input); + Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL; + if(dim==3){ + sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input); + sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input); + sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input); + } + + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsTensor(tbasis,gauss); + + /*Get tau from inputs*/ + sigmapxx_input->GetInputValue(&sigmapxx,gauss); + sigmapyy_input->GetInputValue(&sigmapyy,gauss); + sigmapxy_input->GetInputValue(&sigmapxy,gauss); + if(dim==3){ + sigmapzz_input->GetInputValue(&sigmapzz,gauss); + sigmapxz_input->GetInputValue(&sigmapxz,gauss); + sigmapyz_input->GetInputValue(&sigmapyz,gauss); + } + + /*Calculate d from previous results*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==3){ + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + } + epsxx = dvx[0]; + epsyy = dvy[1]; + epsxy = 0.5*(dvx[1] + dvy[0]); + if(dim==3){ + epszz = dvz[2]; + epsxz = 0.5*(dvx[2] + dvz[0]); + epsyz = 0.5*(dvy[2] + dvz[1]); + } + + /*Solve 2 eta_0 |d|^s-1 + r |d| = |rD(u) + tau|*/ + IssmDouble coef1,coef2,coef3; + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + coef1 = (B/2.)*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 + coef2 = r; + if(dim==2){ + coef3 = sqrt( + (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx) + + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy) + + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy) + ); + } + else{ + coef3 = sqrt( + (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx) + + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy) + + (r*epszz + sigmapzz)*(r*epszz + sigmapzz) + + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy) + + 2*(r*epsxz + sigmapxz)*(r*epsxz + sigmapxz) + + 2*(r*epsyz + sigmapyz)*(r*epsyz + sigmapyz) + ); + } + IssmDouble dnorm; + NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n); + + /*Create Ke*/ + D_scalar=(coef1*pow(dnorm,(1.-n)/n)+r)*gauss->weight*Jdet; + TripleMultiply(tbasis,tnumnodes,1,0, + &D_scalar,1,1,0, + tbasis,1,tnumnodes,0, + Ke,1); + + /*Create Right hand sides*/ + for(int ii=0;iiweight*Jdet; + for(int ii=0;iiweight*Jdet; + for(int ii=0;iiweight*Jdet; + if(dim==3){ + for(int ii=0;iiweight*Jdet; + for(int ii=0;iiweight*Jdet; + for(int ii=0;iiweight*Jdet; + } + } + + /*Solve the systems*/ + _error_("S"); + + /*Clean up and */ + xDelete(tbasis); + xDelete(Ke); + xDelete(pe_xx); + xDelete(pe_yy); + xDelete(pe_zz); + xDelete(pe_xy); + xDelete(pe_xz); + xDelete(pe_yz); } Index: /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp (revision 17546) +++ /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp (revision 17546) @@ -0,0 +1,38 @@ +#include +#include "./types.h" + +int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n){ + /* solve the following equation using Newton-Raphson + * + * c1*x^(s-1) + c2*x = c3 + * + * s = (1+n)/n + * + * we solve y = 10^x: + */ + + + /*Intermediaries*/ + IssmDouble s = (1.+n)/n; + IssmDouble y2; + IssmDouble threshold = 1.e-12; + + /*Initial guess*/ + IssmDouble y1 = 0.; + + while(true){ + + /*Newton step*/ + y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10)*pow(pow(10.,y1),s-1) + c2*log(10.)*pow(10.,y1)); + + if( fabs(y2-y1)/fabs(y1)