Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15424) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15425) @@ -6830,5 +6830,4 @@ switch(fe_stokes){ case 0: - printf("fe %i\n",fe_stokes); /*compute all stiffness matrices for this element*/ Ke1=CreateKMatrixDiagnosticStokesViscous(); @@ -6838,6 +6837,6 @@ case 1: /*compute all stiffness matrices for this element*/ - Ke1=CreateKMatrixDiagnosticStokesViscous(); - Ke2=CreateKMatrixDiagnosticStokesFriction(); + Ke1=CreateKMatrixDiagnosticStokesGLSViscous(); + Ke2=CreateKMatrixDiagnosticStokesGLSFriction(); Ke =new ElementMatrix(Ke1,Ke2); break; @@ -6914,6 +6913,136 @@ } /*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSViscous {{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){ + + /*Intermediaries */ + int i,approximation; + IssmDouble Jdet,viscosity,stokesreconditioning; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble B[8][24]; + IssmDouble B_prime[8][24]; + IssmDouble D_scalar; + IssmDouble D[8][8]={0.0}; + IssmDouble Ke_temp[24][24]={1.0}; //for the six nodes + GaussPenta *gauss=NULL; + + /*If on water or not Stokes, skip stiffness: */ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBStokesGLS(&B[0][0],&xyz_list[0][0],gauss); + GetBprimeStokesGLS(&B_prime[0][0],&xyz_list[0][0],gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + D_scalar=gauss->weight*Jdet; + for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; + for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning; + + TripleMultiply( &B[0][0],8,24,1, + &D[0][0],8,8,0, + &B_prime[0][0],8,24,0, + &Ke_temp[0][0],1); + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/ ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + const int numdof2d=NUMVERTICES2D*NDOF4; + + /*Intermediaries */ + int i,j; + int analysis_type,approximation; + IssmDouble alpha2,Jdet2d; + IssmDouble stokesreconditioning,viscosity; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble LStokes[2][numdof2d]; + IssmDouble DLStokes[2][2]={0.0}; + IssmDouble Ke_drag_gaussian[numdof2d][numdof2d]; + Friction* friction=NULL; + GaussPenta *gauss=NULL; + + /*If on water or not Stokes, skip stiffness: */ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); + GetLStokes(&LStokes[0][0], gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); + + DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx + DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy + + TripleMultiply( &LStokes[0][0],2,numdof2d,1, + &DLStokes[0][0],2,2,0, + &LStokes[0][0],2,numdof2d,0, + &Ke_drag_gaussian[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_drag_gaussian[i][j]; + } + + /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ + //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + + /*Clean up and return*/ + delete gauss; + delete friction; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSFriction{{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSFriction(void){ /*Constants*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15424) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15425) @@ -249,5 +249,7 @@ ElementMatrix* CreateKMatrixDiagnosticStokes(void); ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void); + ElementMatrix* CreateKMatrixDiagnosticStokesGLSViscous(void); ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void); + ElementMatrix* CreateKMatrixDiagnosticStokesGLSFriction(void); ElementMatrix* CreateKMatrixDiagnosticVert(void); ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15424) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15425) @@ -314,4 +314,71 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetBStokesGLS {{{*/ +void PentaRef::GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){ + + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. + * For node i, Bi can be expressed in the actual coordinate system + * by: Bi=[ dh/dx 0 0 0 ] + * [ 0 dh/dy 0 0 ] + * [ 0 0 dh/dy 0 ] + * [ 1/2*dh/dy 1/2*dh/dx 0 0 ] + * [ 1/2*dh/dz 0 1/2*dh/dx 0 ] + * [ 0 1/2*dh/dz 1/2*dh/dy 0 ] + * [ 0 0 0 h ] + * [ dh/dx dh/dy dh/dz 0 ] + * where h is the interpolation function for node i. + * Same thing for Bb except the last column that does not exist. + */ + + int i; + + IssmDouble dh1dh6[3][NUMNODESP1]; + IssmDouble l1l6[NUMNODESP1]; + + /*Get dh1dh7 in actual coordinate system: */ + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); + GetNodalFunctionsP1(l1l6, gauss); + + /*Build B: */ + for (i=0;i