| 118 | } |
| 119 | |
| 120 | int |
| 121 | PFEMElement3Dmini::update() |
| 122 | { |
| 123 | int numnodes = ntags.Size()/2; |
| 124 | int ndm = body.Size(); |
| 125 | |
| 126 | // get nodal coordinates |
| 127 | Vector x(numnodes), y(numnodes), z(numnodes); |
| 128 | for(int i=0; i<numnodes; i++) { |
| 129 | const Vector& coord = nodes[2*i]->getCrds(); |
| 130 | const Vector& disp = nodes[2*i]->getTrialDisp(); |
| 131 | x(i) = coord(0)+disp(0); |
| 132 | y(i) = coord(1)+disp(1); |
| 133 | z(i) = coord(2)+disp(2); |
| 134 | } |
| 135 | |
| 136 | // get Jacobi |
| 137 | J = det(x(1)-x(0),y(1)-y(0),z(1)-z(0), |
| 138 | x(2)-x(0),y(2)-y(0),z(2)-z(0), |
| 139 | x(3)-x(0),y(3)-y(0),z(3)-z(0)); |
| 140 | |
| 141 | if(checkJ && J<=0) { |
| 142 | opserr<<"WARNING: element "<<this->getTag()<<" Jacobian determinant "; |
| 143 | opserr<<J<<" <= 0\n"; |
| 144 | opserr<<"J = "<<J<<"\n"; |
| 145 | opserr<<"Jn = "<<Jn<<"\n"; |
| 146 | opserr<<"x = "<<x; |
| 147 | opserr<<", y = "<<y; |
| 148 | opserr<<", z = "<<z<<"\n"; |
| 149 | return -1; |
| 150 | } |
| 151 | |
| 152 | // deformation tensor |
| 153 | Matrix F1(ndm,ndm), F(ndm,ndm); |
| 154 | dNdx.resize(ndm,numnodes); dNdx.Zero(); |
| 155 | for(int i=0; i<ndm; i++) { |
| 156 | F1(i,0) = x(i+1)-x(0); |
| 157 | F1(i,1) = y(i+1)-y(0); |
| 158 | F1(i,2) = z(i+1)-z(0); |
| 159 | } |
| 160 | F1.Invert(F); |
| 161 | |
| 162 | // get dNdx |
| 163 | for(int i=0; i<ndm; i++) { |
| 164 | for(int a=1; a<numnodes; a++) { |
| 165 | dNdx(i,a) = F(i,a-1)*J; |
| 166 | dNdx(i,0) -= F(i,a-1)*J; |
| 167 | } |
| 168 | } |
| 169 | |
| 170 | #ifdef debug3dmini |
| 171 | Vector xn(numnodes), yn(numnodes), zn(numnodes); |
| 172 | for(int i=0; i<numnodes; i++) { |
| 173 | const Vector& coord = nodes[2*i]->getCrds(); |
| 174 | const Vector& dispn = nodes[2*i]->getDisp(); |
| 175 | xn(i) = coord(0)+dispn(0); |
| 176 | yn(i) = coord(1)+dispn(1); |
| 177 | zn(i) = coord(2)+dispn(2); |