| 922 | } |
| 923 | |
| 924 | void |
| 925 | SSPquad::GetStab(void) |
| 926 | // this function computes the stabilization stiffness matrix for the element |
| 927 | { |
| 928 | Vector g1(SSPQ_NUM_DIM); |
| 929 | Vector g2(SSPQ_NUM_DIM); |
| 930 | Matrix I(SSPQ_NUM_DIM,SSPQ_NUM_DIM); |
| 931 | Matrix FCF(SSPQ_NUM_DIM,SSPQ_NUM_DIM); |
| 932 | Matrix Jmat(SSPQ_NUM_DIM,SSPQ_NUM_DIM); |
| 933 | Matrix Jinv(SSPQ_NUM_DIM,SSPQ_NUM_DIM); |
| 934 | Matrix dNloc(SSPQ_NUM_NODE,SSPQ_NUM_DIM); |
| 935 | Matrix dN(SSPQ_NUM_NODE,SSPQ_NUM_DIM); |
| 936 | Matrix Mben(2,SSPQ_NUM_DOF); |
| 937 | double Hss; |
| 938 | double Hst; |
| 939 | double Htt; |
| 940 | |
| 941 | // shape function derivatives (local crd) at center |
| 942 | dNloc(0,0) = -0.25; |
| 943 | dNloc(1,0) = 0.25; |
| 944 | dNloc(2,0) = 0.25; |
| 945 | dNloc(3,0) = -0.25; |
| 946 | dNloc(0,1) = -0.25; |
| 947 | dNloc(1,1) = -0.25; |
| 948 | dNloc(2,1) = 0.25; |
| 949 | dNloc(3,1) = 0.25; |
| 950 | |
| 951 | // jacobian matrix |
| 952 | Jmat = mNodeCrd*dNloc; |
| 953 | // inverse of the jacobian matrix |
| 954 | Jmat.Invert(Jinv); |
| 955 | |
| 956 | // shape function derivatives (global crd) |
| 957 | dN = dNloc*Jinv; |
| 958 | |
| 959 | // define hourglass stabilization vector gamma = 0.25*(h - (h^x)*bx - (h^y)*by); |
| 960 | double hx = mNodeCrd(0,0) - mNodeCrd(0,1) + mNodeCrd(0,2) - mNodeCrd(0,3); |
| 961 | double hy = mNodeCrd(1,0) - mNodeCrd(1,1) + mNodeCrd(1,2) - mNodeCrd(1,3); |
| 962 | double gamma[4]; |
| 963 | gamma[0] = 0.25*( 1.0 - hx*dN(0,0) - hy*dN(0,1)); |
| 964 | gamma[1] = 0.25*(-1.0 - hx*dN(1,0) - hy*dN(1,1)); |
| 965 | gamma[2] = 0.25*( 1.0 - hx*dN(2,0) - hy*dN(2,1)); |
| 966 | gamma[3] = 0.25*(-1.0 - hx*dN(3,0) - hy*dN(3,1)); |
| 967 | |
| 968 | // define mapping matrices |
| 969 | Mmem.Zero(); |
| 970 | Mben.Zero(); |
| 971 | for (int i = 0; i < 4; i++) { |
| 972 | Mmem(0,2*i) = dN(i,0); |
| 973 | Mmem(1,2*i+1) = dN(i,1); |
| 974 | Mmem(2,2*i) = dN(i,1); |
| 975 | Mmem(2,2*i+1) = dN(i,0); |
| 976 | |
| 977 | Mben(0,2*i) = gamma[i]; |
| 978 | Mben(1,2*i+1) = gamma[i]; |
| 979 | } |
| 980 | |
| 981 | // base vectors |
nothing calls this directly
no test coverage detected