| 158 | } |
| 159 | |
| 160 | int |
| 161 | PFEMElement2D::update() |
| 162 | { |
| 163 | // get nodal coordinates |
| 164 | double x[3], y[3]; |
| 165 | for(int a=0; a<3; a++) { |
| 166 | const Vector& coord = nodes[2*a]->getCrds(); |
| 167 | const Vector& disp = nodes[2*a]->getTrialDisp(); |
| 168 | x[a] = coord(0) + disp(0); |
| 169 | y[a] = coord(1) + disp(1); |
| 170 | } |
| 171 | |
| 172 | // get c and d |
| 173 | double cc[3], dd[3]; |
| 174 | cc[0] = y[1]-y[2]; |
| 175 | dd[0] = x[2]-x[1]; |
| 176 | cc[1] = y[2]-y[0]; |
| 177 | dd[1] = x[0]-x[2]; |
| 178 | cc[2] = y[0]-y[1]; |
| 179 | dd[2] = x[1]-x[0]; |
| 180 | |
| 181 | // get Jacobi |
| 182 | double J = cc[0]*dd[1]-dd[0]*cc[1]; |
| 183 | |
| 184 | if(fabs(J)<1e-15) { |
| 185 | opserr<<"WARNING: element area is nearly zero"; |
| 186 | opserr<<" -- PFEMElement2D::update\n"; |
| 187 | for (int i=0; i<3; i++) { |
| 188 | opserr<<"node "<<ntags[2*i]<<"\n"; |
| 189 | opserr<<"x = "<<x[i]<<" , y = "<<y[i]<<"\n"; |
| 190 | } |
| 191 | |
| 192 | return -1; |
| 193 | } |
| 194 | |
| 195 | // get M |
| 196 | M = rho*J*thickness/6.0; |
| 197 | Mp = (kappa<=0? 0.0 : J*thickness/kappa/24.0); |
| 198 | double Mb = 9.*rho*J*thickness/40.0; |
| 199 | |
| 200 | // get Km |
| 201 | Km.Zero(); |
| 202 | double fact = mu*thickness/(6.*J); |
| 203 | for (int a=0; a<3; a++) { |
| 204 | for (int b=0; b<3; b++) { |
| 205 | Km(2*a,2*b) = fact*(4*cc[a]*cc[b]+3*dd[a]*dd[b]); // Kxx |
| 206 | Km(2*a,2*b+1) = fact*(3*dd[a]*cc[b]-2*cc[a]*dd[b]); // Kxy |
| 207 | Km(2*a+1,2*b) = fact*(3*cc[a]*dd[b]-2*dd[a]*cc[b]); // Kyx |
| 208 | Km(2*a+1,2*b+1) = fact*(3*cc[a]*cc[b]+4*dd[a]*dd[b]); // Kyy |
| 209 | } |
| 210 | } |
| 211 | |
| 212 | // get Kb |
| 213 | Matrix Kb(2,2); |
| 214 | fact = 27.*mu*thickness/(40.*J); |
| 215 | double cc2 = 0., dd2 = 0., cd2 = 0.; |
| 216 | for(int a=0; a<3; a++) { |
| 217 | cc2 += cc[a]*cc[a]; |
nothing calls this directly
no test coverage detected