1 /* 2 * This program is free software; you can redistribute it and/or 3 * modify it under the terms of the GNU General Public License as 4 * published by the Free Software Foundation; either version 2 of the 5 * License, or (at your option) any later version. 6 * 7 * This program is distributed in the hope that it will be useful, but 8 * WITHOUT ANY WARRANTY; without even the implied warranty of 9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOUSE. See the GNU 10 * General Public License for more details. 11 * 12 * You should have recieved a copy of the GNU General Public License 13 * along with this program; if not write to the Free Software 14 * Foundation, inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 15 * USA 16 */ 17 18 package run.elements; 19 20 import run.*; 21 22 // FILE: c:/temp/krockpackage/Solid.java 23 import java.util.*; 24 25 26 /** 27 * A class that represents ... 28 * 29 * @author Jonas Forssell, Yuriy Mikhaylovskiy 30 * 31 * @see OtherClasses 32 */ 33 public class Solid_Iso_6HG extends Element { 34 private Material[] material; 35 private Jama.Matrix D; 36 private int number_of_integration_points; 37 private double xsi; 38 private double phi; 39 private double etha; 40 private double AHR; 41 private Jama.Matrix H; 42 private Jama.Matrix M; 43 private Jama.Matrix N; 44 private Jama.Matrix B; 45 private Jama.Matrix J; 46 private Jama.Matrix J_inv; 47 private Jama.Matrix strain; 48 private Jama.Matrix dstrain; 49 private Jama.Matrix stress; 50 private Jama.Matrix[] hforce; 51 52 private double [][] elemp; 53 54 private Jama.Matrix d; 55 private Jama.Matrix f; 56 private Jama.Matrix P; 57 private double[] W; 58 private boolean AHR_is_set; 59 private boolean NIP_is_set; 60 private boolean Nodes_are_set; 61 private boolean Material_is_set; 62 63 /** 64 * This is the Solid Isoparametric element. The element can have one or 65 * eight integration points. Creation date: (2001-08-10 19:46:30) 66 */ Solid_Iso_6HG()67 public Solid_Iso_6HG() { 68 super(); 69 type = new String("SOLID_ISO_6HG"); 70 71 int i; 72 73 // 74 material = new Material[1]; 75 76 // 77 xsi = 0; 78 etha = 0; 79 phi = 0; 80 node = new Node[8]; 81 W = new double[1]; 82 H = new Jama.Matrix(6, 9); 83 M = new Jama.Matrix(8, 3); 84 d = new Jama.Matrix(24, 1); 85 f = new Jama.Matrix(24, 1); 86 87 B = new Jama.Matrix(6, 24); 88 89 D = new Jama.Matrix(3, 8); 90 91 J = new Jama.Matrix(3, 3); 92 93 J_inv = new Jama.Matrix(3, 3); 94 95 P = new Jama.Matrix(9, 9); 96 N = new Jama.Matrix(9, 24); 97 98 strain = new Jama.Matrix(6, 1); 99 // 100 dstrain = new Jama.Matrix(6, 1); 101 // 102 stress = new Jama.Matrix(6, 1); 103 104 hforce = new Jama.Matrix[8]; 105 106 for (i = 0; i < 8; i++) { 107 hforce[i] = new Jama.Matrix(3, 1); 108 } 109 110 111 112 // Setting default value for number of integration points 113 number_of_integration_points = 1; 114 } 115 116 /** 117 * This method calculates the mass matrix of the element. The matrix is 118 * lumped, i.e. the mass is concentrated to the element nodes. There is a 119 * consistent way of determining the mass distribution by using the HRZ 120 * Lumping Scheme 1. Compute the diagonal coefficients of the consistent 121 * mass matrix 2. Compute the total mass of the element m 3. Compute a 122 * number s by adding the diagonal coefficients Mii associated with 123 * translational d.o.f (but not rotational d.o.f, if any) that are 124 * mutually parallell and in the same direction 4. Scale all the diagonal 125 * coefficients by multiplying them by the ratio m/s, thus preserving the 126 * total mass of the element. The mass matrix m = integral 127 * (densityN_transposeNdV) This integral can be translated into a gauss 128 * point summation. Creation date: (25/12/01 %T) 129 */ assembleMassMatrix()130 public void assembleMassMatrix() 131 throws IllegalArgumentException 132 { 133 int i; 134 int j; 135 int k; 136 Jama.Matrix mass; 137 double total_mass; 138 double s; 139 double det_j; 140 141 // Initialize mass 142 mass = new Jama.Matrix(24, 24); 143 total_mass = 0; 144 s = 0; 145 146 /* This element has 8 nodes and with 3 direction each. This gives a 24x24 matrix. 147 * The matrix will be different for each integration point. */ 148 149 // Define M (the node coordinate matrix) 150 for (j = 0; j < 8; j++) { 151 M.set(j, 0, node[j].getX_pos()); 152 M.set(j, 1, node[j].getY_pos()); 153 M.set(j, 2, node[j].getZ_pos()); 154 } 155 156 // Now, start by computing N, via D 157 // Compute D. Use only one integration point. Then xsi=phi=etha=0 158 calculateD(0, 0, 0); 159 160 // Create the N matrix by expanding the D matrix 161 calculateN(); 162 163 // Now, do the matrix calculation to derive the jacobian. 164 J = D.times(M); 165 det_j = J.det(); 166 167 // Sum up the mass contribution (The weight function for one gauss point is 2^3 = 8) 168 mass = N.transpose().times(N).times(material[0].getDensity()) 169 .times(det_j).times(8.0); 170 171 // Keep only the diagonal elements 172 for (j = 0; j < 24; j++) { 173 for (k = 0; k < 24; k++) { 174 mass.set(j, k, (j == k) ? mass.get(j, k) : 0); 175 } 176 } 177 178 179 /* Now, we have the 24x24 mass matrix. Move on to step 2, compute total mass of the element. 180 * The total mass is equal to the volume integrand of density. Gauss translation means 181 * sum[density*jacobian*weightfunction] over all integration points.*/ 182 // Compute D 183 //calculateD(0, 0, 0); 184 185 // Now, do the matrix calculation to derive the jacobian. 186 //J = D.times(M); 187 188 // Sum up the mass contribution (Weight function 2^3 = 8) 189 total_mass = material[0].getDensity() * det_j * 8; 190 191 // Check mass for element 192 if (total_mass <= 0) { 193 throw new IllegalArgumentException( 194 "Error in Solid Hexahedron Element " + number + 195 ". Element mass is zero or negative.\nMaterial density: " + 196 material[0].getDensity() 197 ); 198 } 199 200 if (det_j <= 0) { 201 throw new IllegalArgumentException( 202 "Error in Solid Hexahedron Element " + number + 203 ". Element Volume is zero or negative. Check node defintion in indata file." 204 ); 205 } 206 207 // Third stage is to compute a scaling factor s using u-direction for all nodes 208 for (i = 0; i < 8; i++) { 209 s += mass.get(3 * i, 0); 210 } 211 212 //Fourth stage is to scale all diagonal coefficients using total_mass/s as a scaling factor 213 for (i = 0; i < 24; i++) { 214 mass.set(i, i, (mass.get(i, i) * total_mass) / s); 215 } 216 217 218 // We now have a lumped mass matrix stored in the mass matrix. It must be distributed to the nodes. 219 // Use only the u-component. Assume the others are the same. (They should be) 220 for (i = 0; i < 8; i++) { 221 node[i].addMass(mass.get(3 * i, 3 * i) / 8.0); 222 } 223 224 // That's it! Finish. 225 } 226 227 /** 228 * Insert the method's description here. Creation date: (25/12/01 %T) 229 */ calculateContactForces()230 public void calculateContactForces() { 231 } 232 233 /** 234 * Calculates the D-matrix for integration point i Creation date: (05/01/02 235 * %T) 236 * 237 * @param integration_point int 238 */ calculateD(double xsi, double phi, double etha)239 private void calculateD(double xsi, double phi, double etha) { 240 // Calculate the D-matrix for integration point i 241 // WAS CHANGED TO REFLECT PARAMETRIC COORDINATES OF HOURGLASS 242 D.set(0, 0, -(1.0 / 8) ); 243 D.set(0, 1, +(1.0 / 8) ); 244 D.set(0, 2, +(1.0 / 8) ); 245 D.set(0, 3, -(1.0 / 8) ); 246 D.set(0, 4, -(1.0 / 8) ); 247 D.set(0, 5, +(1.0 / 8) ); 248 D.set(0, 6, +(1.0 / 8) ); 249 D.set(0, 7, -(1.0 / 8) ); 250 251 // 252 D.set(1, 0, -(1.0 / 8) ); 253 D.set(1, 1, -(1.0 / 8) ); 254 D.set(1, 2, (1.0 / 8) ); 255 D.set(1, 3, (1.0 / 8) ); 256 D.set(1, 4, -(1.0 / 8) ); 257 D.set(1, 5, -(1.0 / 8) ); 258 D.set(1, 6, (1.0 / 8) ); 259 D.set(1, 7, (1.0 / 8) ); 260 261 // 262 D.set(2, 0, -(1.0 / 8) ); 263 D.set(2, 1, -(1.0 / 8) ); 264 D.set(2, 2, -(1.0 / 8) ); 265 D.set(2, 3, -(1.0 / 8) ); 266 D.set(2, 4, (1.0 / 8) ); 267 D.set(2, 5, (1.0 / 8) ); 268 D.set(2, 6, (1.0 / 8) ); 269 D.set(2, 7, (1.0 / 8) ); 270 } 271 272 /** 273 * Insert the method's description here. Creation date: (25/12/01 %T) 274 */ calculateExternalForces(double currtime)275 public void calculateExternalForces(double currtime) { 276 } 277 278 /** 279 * Calculates the N-matrix for integration point i Creation date: (05/01/02 280 * %T) 281 */ calculateN()282 private void calculateN() { 283 // Calculate the N-matrix for integration point i by expanding the D-matrix 284 N.setMatrix(0, 2, 0, 0, D.getMatrix(0, 2, 0, 0)); 285 N.setMatrix(3, 5, 1, 1, D.getMatrix(0, 2, 0, 0)); 286 N.setMatrix(6, 8, 2, 2, D.getMatrix(0, 2, 0, 0)); 287 N.setMatrix(0, 2, 3, 3, D.getMatrix(0, 2, 1, 1)); 288 N.setMatrix(3, 5, 4, 4, D.getMatrix(0, 2, 1, 1)); 289 N.setMatrix(6, 8, 5, 5, D.getMatrix(0, 2, 1, 1)); 290 N.setMatrix(0, 2, 6, 6, D.getMatrix(0, 2, 2, 2)); 291 N.setMatrix(3, 5, 7, 7, D.getMatrix(0, 2, 2, 2)); 292 N.setMatrix(6, 8, 8, 8, D.getMatrix(0, 2, 2, 2)); 293 N.setMatrix(0, 2, 9, 9, D.getMatrix(0, 2, 3, 3)); 294 N.setMatrix(3, 5, 10, 10, D.getMatrix(0, 2, 3, 3)); 295 N.setMatrix(6, 8, 11, 11, D.getMatrix(0, 2, 3, 3)); 296 N.setMatrix(0, 2, 12, 12, D.getMatrix(0, 2, 4, 4)); 297 N.setMatrix(3, 5, 13, 13, D.getMatrix(0, 2, 4, 4)); 298 N.setMatrix(6, 8, 14, 14, D.getMatrix(0, 2, 4, 4)); 299 N.setMatrix(0, 2, 15, 15, D.getMatrix(0, 2, 5, 5)); 300 N.setMatrix(3, 5, 16, 16, D.getMatrix(0, 2, 5, 5)); 301 N.setMatrix(6, 8, 17, 17, D.getMatrix(0, 2, 5, 5)); 302 N.setMatrix(0, 2, 18, 18, D.getMatrix(0, 2, 6, 6)); 303 N.setMatrix(3, 5, 19, 19, D.getMatrix(0, 2, 6, 6)); 304 N.setMatrix(6, 8, 20, 20, D.getMatrix(0, 2, 6, 6)); 305 N.setMatrix(0, 2, 21, 21, D.getMatrix(0, 2, 7, 7)); 306 N.setMatrix(3, 5, 22, 22, D.getMatrix(0, 2, 7, 7)); 307 N.setMatrix(6, 8, 23, 23, D.getMatrix(0, 2, 7, 7)); 308 } 309 310 /** 311 * This method will calculate the nodal forces resulting from the stresses 312 * in the element. The forces are all added up in a force matrix f with 313 * the following format: f = [ fx1 fy1 fz1 fx2 fy2 ..... fz8 ] Creation 314 * date: (25/12/01 %T) 315 */ calculateNodalForces(int dummy, double dt)316 public void calculateNodalForces(int dummy, double dt) { 317 int n; 318 int i, j, k, l; 319 Jama.Matrix global_force; 320 double A = 0; 321 double vol = 0; 322 double S0 = 0; 323 double G = 0; 324 double X[][]; 325 double V[][]; 326 double BB[][]; 327 double GS[][]= {{ 1., 1.,-1.,-1.,-1.,-1., 1., 1.}, 328 { 1.,-1.,-1., 1.,-1., 1., 1.,-1.}, 329 { 1.,-1., 1.,-1., 1.,-1., 1.,-1.}, 330 {-1., 1.,-1., 1., 1.,-1., 1.,-1.}}; 331 double GB[][]= {{ 1., 1.,-1.,-1.,-1.,-1., 1., 1.}, 332 { 1.,-1.,-1., 1.,-1., 1., 1.,-1.}, 333 { 1.,-1., 1.,-1., 1.,-1., 1.,-1.}, 334 {-1., 1.,-1., 1., 1.,-1., 1.,-1.}}; 335 336 337 /* The nodal force array is the sum of all integration point contributions times the gauss weight factor. 338 * In this particular element, where the order of gauss evaluation is 2, the weight factor is 1 and the 339 * multiplication is therefore not done. Otherwise, the expression would have been: 340 * 341 * f = B.transpose().times(stress[i]).times(J.det()).times(weightfactor); 342 * 343 * Since this method is called once for every integration point, this will happen. 344 * Start by calculating the force vector contribution */ 345 f = B.transpose().times(stress).times(J.det()).times(8); 346 347 /* Now, all the node contributions are represented in this matrix. It needs to be split up. 348 * Each node contribution must then be added to each node 349 */ 350 351 352 for (n = 0; n < 8; n++) { 353 // Split up; extract a 3x1 matrix from the large matrix. global_force = [ fxN fyN fzN ] 354 global_force = f.getMatrix(3 * n, (3 * n) + 2, 0, 0); 355 356 // Add to node (Internal forces are subtracted) 357 node[n].addInternalForce(global_force.times(-1.0)); 358 } 359 360 // the jama is realy time consuming and ineficient for small mtrices... 361 // I prefer arrays, instead... Sorry. 362 BB = B.getArray(); 363 vol = 8*J.det(); 364 365 /* 366 367 GS=GB; // I don't know if this will work, but... 368 for(i=0; i<8; i++) { 369 for(j=0; i<3; j++) { 370 A=A+BB[i][j]*BB[i][j]; 371 S0=BB[i][j]/vol; 372 for(k=0; k<4; k++) { 373 //GS[i][k]=GB[i][k]; 374 for(l=0; l<8; l++) { 375 GS[i][k] = GS[i][k] - S0 * X[l][j]*GB[l][k]; 376 } 377 } 378 } 379 } 380 381 382 // Shear modulus... Need to be used here... Please help in doing so... 383 G= (material[0].getYoungsModulus())/(2*(1+(material[0].getNu()))); 384 385 // Global time step, dt, is needed here... How to get it? 386 // AHR is an user defined parameter to control HG... we need it together 387 // with the material properties... 388 A= dt*AHR*(material[0].getYoungsModulus()+G+G)*A/(24*vol); 389 for(i=0; i<3; i++) { 390 for(k=0; k<4; k++) { 391 S0=0.; 392 for( j=0; j<8; j++) { 393 S0= S0 + GS[j][k]*V[i][j]; 394 } 395 S0=elemp[i][k]+A*S0; 396 elemp[i][k]=S0; 397 for( j=0; j<8; j++) { 398 // We need to add the hour glass stabilization force to 399 // the global force vector... The indices are a little confusing 400 // but it works... 401 // i is the direction (x, y or z) 402 // j is the local node number 403 // k is the assiciated HG mode 404 // F[i][j]=F[i][j]+S0*GS[j][k]; 405 hforce[j].set(i,1,hforce[j].get(i,1)+S0*GS[j][k]); 406 } 407 } 408 } 409 410 for (i=0; i<8; i++) 411 node[i].addHourglassForce(hforce[i]); 412 */ 413 } 414 415 /** 416 * Definition of strain vector is [exx eyy ezz gammaxy gammayz gammxz] The 417 * strain is derived as strain = B d where d is an array of node 418 * displacements d = [u1 v1 w1 u2 v2 w2 .... w8] The B-matrix is derived 419 * from three main matrices, H, J and D according to: strain = HPN where P 420 * expanded version of the inverted jacobian J N is an expanded version of 421 * D J is called the jacobian. It consists of an array M which includes 422 * the node coordinates and the D matrix again as: J = DM H is a coupling 423 * matrix. H is defined in the setInitialConditions method since it is not 424 * dependent on any variable. Note, this method is run several times. i is 425 * the number of the integration point currently run. Creation date: 426 * (25/12/01 %T) 427 * 428 * @param tstep double 429 */ calculateStrain(double tstep, int i)430 public void calculateStrain(double tstep, int i) { 431 int j; 432 433 // Calculate the D-matrix 434 calculateD(xsi, phi, etha); 435 436 // Create the N matrix by expanding the D matrix 437 calculateN(); 438 439 // Update M (the node coordinate matrix) 440 for (j = 0; j < 8; j++) { 441 M.set(j, 0, node[j].getX_pos()); 442 M.set(j, 1, node[j].getY_pos()); 443 M.set(j, 2, node[j].getZ_pos()); 444 } 445 446 // And determine the displacement matrix (difference between old, one timestep ago and new node positions) 447 for (j = 0; j < 8; j++) { 448 d.set(3 * j, 0, node[j].getX_dpos()); 449 d.set((3 * j) + 1, 0, node[j].getY_dpos()); 450 d.set((3 * j) + 2, 0, node[j].getZ_dpos()); 451 } 452 453 // Now, do the matrix calculation to derive the jacobian. 454 J = D.times(M); 455 456 // Now, create the P matrix by first inverting it.. 457 J_inv = J.inverse(); 458 459 // ... and then expanding it 460 P.setMatrix(0, 2, 0, 2, J_inv); 461 P.setMatrix(3, 5, 3, 5, J_inv); 462 P.setMatrix(6, 8, 6, 8, J_inv); 463 464 // Now, calculate the B matrix 465 B = H.times(P.times(N)); 466 467 // Finally, calculate the strain increment array for this integration point 468 dstrain = B.times(d); 469 } 470 471 /** 472 * This method calculates the stresses in each integration point based on 473 * old strain, strain increment and old stress. The material law will 474 * update the strain and stress matrices to new values using constitutive 475 * law and strain increment dstrain. The stresses are: [Sxx Syy Szz Sxy 476 * Syz Sxz]T (or [sigmaxx sigmayy sigmazz tauxy tauyz tauxz]T). Note tauxz 477 * = tauzx and so on. Creation date: (25/12/01 %T) 478 */ calculateStress(int i, double timestep)479 public void calculateStress(int i, double timestep) { 480 // i is the current integration point. 481 material[i].calculateStressThreeDimensional( 482 strain, dstrain, stress, timestep 483 ); 484 } 485 486 /** 487 * The critical timestep for a solid element is related to the minimum 488 * distance through the element. This can normally be found by dividing 489 * the element volume by the largest side area. However, in this special 490 * case with a hexahedron, the easiest way is to to calculate the shortest 491 * edge length since this is always the critical length. Creation date: 492 * (25/12/01 %T) 493 * 494 * @param current_timestep double 495 * 496 * @return double 497 */ checkTimestep(double current_timestep)498 public double checkTimestep(double current_timestep) { 499 double critical_length; 500 double timestep; 501 double c; 502 503 // Lower edges 504 critical_length = Math.sqrt( 505 ((node[1].getX_pos() - node[0].getX_pos()) * (node[1].getX_pos() - 506 node[0].getX_pos())) + 507 ((node[1].getY_pos() - node[0].getY_pos()) * (node[1].getY_pos() - 508 node[0].getY_pos())) + 509 ((node[1].getZ_pos() - node[0].getZ_pos()) * (node[1].getZ_pos() - 510 node[0].getZ_pos())) 511 ); 512 513 c = Math.sqrt( 514 ((node[2].getX_pos() - node[1].getX_pos()) * (node[2].getX_pos() - 515 node[1].getX_pos())) + 516 ((node[2].getY_pos() - node[1].getY_pos()) * (node[2].getY_pos() - 517 node[1].getY_pos())) + 518 ((node[2].getZ_pos() - node[1].getZ_pos()) * (node[2].getZ_pos() - 519 node[1].getZ_pos())) 520 ); 521 522 if (c > 0) critical_length = Math.min( 523 critical_length,c 524 ); 525 526 527 c = Math.sqrt( 528 ((node[3].getX_pos() - node[2].getX_pos()) * (node[3].getX_pos() - 529 node[2].getX_pos())) + 530 ((node[3].getY_pos() - node[2].getY_pos()) * (node[3].getY_pos() - 531 node[2].getY_pos())) + 532 ((node[3].getZ_pos() - node[2].getZ_pos()) * (node[3].getZ_pos() - 533 node[2].getZ_pos())) 534 ); 535 536 if (c > 0) critical_length = Math.min( 537 critical_length,c 538 ); 539 540 541 c = Math.sqrt( 542 ((node[0].getX_pos() - node[3].getX_pos()) * (node[0].getX_pos() - 543 node[3].getX_pos())) + 544 ((node[0].getY_pos() - node[3].getY_pos()) * (node[0].getY_pos() - 545 node[3].getY_pos())) + 546 ((node[0].getZ_pos() - node[3].getZ_pos()) * (node[0].getZ_pos() - 547 node[3].getZ_pos())) 548 ); 549 550 // Upper edges 551 if (c > 0) critical_length = Math.min( 552 critical_length,c 553 ); 554 555 c = Math.sqrt( 556 ((node[5].getX_pos() - node[4].getX_pos()) * (node[5].getX_pos() - 557 node[4].getX_pos())) + 558 ((node[5].getY_pos() - node[4].getY_pos()) * (node[5].getY_pos() - 559 node[4].getY_pos())) + 560 ((node[5].getZ_pos() - node[4].getZ_pos()) * (node[5].getZ_pos() - 561 node[4].getZ_pos())) 562 ); 563 564 if (c > 0) critical_length = Math.min( 565 critical_length,c 566 ); 567 568 c = Math.sqrt( 569 ((node[6].getX_pos() - node[5].getX_pos()) * (node[6].getX_pos() - 570 node[5].getX_pos())) + 571 ((node[6].getY_pos() - node[5].getY_pos()) * (node[6].getY_pos() - 572 node[5].getY_pos())) + 573 ((node[6].getZ_pos() - node[5].getZ_pos()) * (node[6].getZ_pos() - 574 node[5].getZ_pos())) 575 ); 576 577 if (c > 0) critical_length = Math.min( 578 critical_length,c 579 ); 580 581 c = Math.sqrt( 582 ((node[7].getX_pos() - node[6].getX_pos()) * (node[7].getX_pos() - 583 node[6].getX_pos())) + 584 ((node[7].getY_pos() - node[6].getY_pos()) * (node[7].getY_pos() - 585 node[6].getY_pos())) + 586 ((node[7].getZ_pos() - node[6].getZ_pos()) * (node[7].getZ_pos() - 587 node[6].getZ_pos())) 588 ); 589 590 if (c > 0) critical_length = Math.min( 591 critical_length,c 592 ); 593 594 c = Math.sqrt( 595 ((node[4].getX_pos() - node[7].getX_pos()) * (node[4].getX_pos() - 596 node[7].getX_pos())) + 597 ((node[4].getY_pos() - node[7].getY_pos()) * (node[4].getY_pos() - 598 node[7].getY_pos())) + 599 ((node[4].getZ_pos() - node[7].getZ_pos()) * (node[4].getZ_pos() - 600 node[7].getZ_pos())) 601 ); 602 603 if (c > 0) critical_length = Math.min( 604 critical_length,c 605 ); 606 607 // Vertical edges 608 c = Math.sqrt( 609 ((node[4].getX_pos() - node[0].getX_pos()) * (node[4].getX_pos() - 610 node[0].getX_pos())) + 611 ((node[4].getY_pos() - node[0].getY_pos()) * (node[4].getY_pos() - 612 node[0].getY_pos())) + 613 ((node[4].getZ_pos() - node[0].getZ_pos()) * (node[4].getZ_pos() - 614 node[0].getZ_pos())) 615 ); 616 617 if (c > 0) critical_length = Math.min( 618 critical_length,c 619 ); 620 621 c = Math.sqrt( 622 ((node[5].getX_pos() - node[1].getX_pos()) * (node[5].getX_pos() - 623 node[1].getX_pos())) + 624 ((node[5].getY_pos() - node[1].getY_pos()) * (node[5].getY_pos() - 625 node[1].getY_pos())) + 626 ((node[5].getZ_pos() - node[1].getZ_pos()) * (node[5].getZ_pos() - 627 node[1].getZ_pos())) 628 ); 629 630 if (c > 0) critical_length = Math.min( 631 critical_length,c 632 ); 633 634 c = Math.sqrt( 635 ((node[6].getX_pos() - node[2].getX_pos()) * (node[6].getX_pos() - 636 node[2].getX_pos())) + 637 ((node[6].getY_pos() - node[2].getY_pos()) * (node[6].getY_pos() - 638 node[2].getY_pos())) + 639 ((node[6].getZ_pos() - node[2].getZ_pos()) * (node[6].getZ_pos() - 640 node[2].getZ_pos())) 641 ); 642 643 if (c > 0) critical_length = Math.min( 644 critical_length,c 645 ); 646 647 c = Math.sqrt( 648 ((node[7].getX_pos() - node[3].getX_pos()) * (node[7].getX_pos() - 649 node[3].getX_pos())) + 650 ((node[7].getY_pos() - node[3].getY_pos()) * (node[7].getY_pos() - 651 node[3].getY_pos())) + 652 ((node[7].getZ_pos() - node[3].getZ_pos()) * (node[7].getZ_pos() - 653 node[3].getZ_pos())) 654 ); 655 656 if (c > 0) critical_length = Math.min( 657 critical_length,c 658 ); 659 660 // 661 timestep = critical_length / material[0].wavespeedThreeDimensional( 662 0.0, 0.0 663 ); 664 665 // Now, return the smallest of the suggested and the calculated timestep. 666 // The smallest element will decide the timestep for the whole simulation. 667 return Math.min(timestep, current_timestep); 668 } 669 670 /** 671 * Insert the method's description here. Creation date: (26/12/01 %T) 672 * 673 * @return int 674 */ getNumberOfIntegrationPoints()675 public int getNumberOfIntegrationPoints() { 676 return number_of_integration_points; 677 } 678 679 /** 680 * Insert the method's description here. Creation date: (25/12/01 %T) 681 * 682 * @param arg1 java.lang.String 683 * @param arg2 java.lang.String 684 * @param arg3 java.lang.String 685 * @param lineno int 686 * @param nodelist java.util.Vector 687 * @param materiallist java.util.Vector 688 */ parse_Fembic( Token[] param, int lineno, RplVector nodelist, RplVector materiallist, RplVector loadlist, Hashtable nodetable )689 public void parse_Fembic( 690 Token[] param, int lineno, RplVector nodelist, RplVector materiallist, 691 RplVector loadlist, Hashtable nodetable 692 ) 693 throws java.text.ParseException 694 { 695 int nodenumber; 696 int j; 697 int i = 0; 698 699 while (i < param.length) { 700 // The nodes of the element are defined 701 if ( 702 param[i].getw().toUpperCase().equals("NODES") && 703 param[i + 1].getw().toUpperCase().equals("=") 704 ) { 705 // Assume now that the nodes are delivered in param3, with the format 706 // [nodenr,nodenr,nodenr,nodenr,nodenr,nodenr,nodenr,nodenr] 707 if ( 708 ! param[i + 2].getw().toUpperCase().startsWith("[") || 709 ! param[i + 2].getw().toUpperCase().endsWith("]") 710 ) { 711 throw new java.text.ParseException( 712 "Error, node number definition should be [nodenr1,nodenr2,nodenr3,nodenr4,nodenr5,nodenr6,nodenr7,nodenr8]", 713 lineno 714 ); 715 } 716 717 // Ok, now find the numbers 718 try { 719 for (j = 0; j < 8; j++) { 720 node[j] = super.findNode( 721 super.getNodeNumber(j + 1, param[i + 2].getw()), 722 nodetable 723 ); 724 } 725 } catch (IllegalArgumentException e) { 726 throw new java.text.ParseException(e.getMessage(), lineno); 727 } 728 729 i += 3; 730 Nodes_are_set = true; 731 } else 732 // The material of the element is defined 733 if ( 734 param[i].getw().toUpperCase().equals("MATERIAL") && 735 param[i + 1].getw().toUpperCase().equals("=") 736 ) { 737 // Assume now that the material name is delivered in param3 738 try { 739 // We want the handle to the material object. 740 material[0] = super.findMaterial( 741 param[i + 2].getw().toUpperCase(), materiallist 742 ); 743 i += 3; 744 Material_is_set = true; 745 } catch (IllegalArgumentException e) { 746 throw e; 747 } 748 } else 749 // The cross section area of the element is defined 750 if ( 751 param[i].getw().toUpperCase().equals("NIP") && 752 param[i + 1].getw().toUpperCase().equals("=") 753 ) { 754 // The value of the cross section area is in param4. Set this in the element 755 number_of_integration_points = (int) param[i + 2].getn(); 756 i += 3; 757 NIP_is_set = true; 758 } else 759 // The cross section area of the element is defined 760 if ( 761 param[i].getw().toUpperCase().equals("AHR") && 762 param[i + 1].getw().toUpperCase().equals("=") 763 ) { 764 // The value of the cross section area is in param4. Set this in the element 765 AHR = param[i + 2].getn(); 766 i += 3; 767 AHR_is_set = true; 768 } else { 769 // Neither material or nodes are defined. Then the parameter is wrong. 770 throw new java.text.ParseException( 771 "Unknown Solid (Hexahedran) element parameter ", lineno 772 ); 773 } 774 } 775 } 776 777 /** 778 * Checks that all mandatory variables have been set 779 */ checkIndata()780 public void checkIndata() 781 throws IllegalArgumentException 782 { 783 // Check that all the required parameters have been parsed 784 if (! Nodes_are_set) { 785 throw new IllegalArgumentException( 786 "No nodes defined for Solid_Iso_6 element nr" + number 787 ); 788 } 789 790 if (! Material_is_set) { 791 throw new IllegalArgumentException( 792 "No Material defined for Solid_Iso_6 element nr" + number 793 ); 794 } 795 796 if (! AHR_is_set) { 797 throw new IllegalArgumentException( 798 "No Hourglass parameter AHR set for Solid_Iso_6HG element nr" + number 799 ); 800 } 801 } 802 803 /** 804 * 805 * The method checks destroyed element or not. 806 * If the element is destroyed variable failed = true else failed = false. 807 * 808 */ checkIfFailed()809 public void checkIfFailed() { 810 } 811 812 /** 813 * This method is used to create the lines needed in the result file. The 814 * method returns a string which is printed directly. However, due to the 815 * fact that the line may be different depending on what is requested to 816 * be printed and that the number of methods should be kept down, the 817 * first parameter here is a control parameter. This parameter describes 818 * what should be printed. The second parameter is a required input when 819 * gauss point results are to be printed. Creation date: (09/12/01 %T) 820 * 821 * @param ctrl int 822 * @param gpn int 823 * 824 * @return java.lang.String 825 */ print_Gid(int ctrl, int gpn)826 public String print_Gid(int ctrl, int gpn) { 827 String out; 828 829 switch (ctrl) { 830 case MESH_HEADER: 831 832 /* Print the header for the mesh. 833 * In this case, the type of element is Hexahedra (6-sided solid element) and it uses 8 nodes 834 */ 835 out = new String( 836 "MESH \"MeshType" + type + 837 "\" Dimension 3 ElemType Hexahedra Nnode 8\n" 838 ); 839 840 return out; 841 842 case MESH: 843 844 /* Print the element number and connected nodes */ 845 out = new String( 846 number + "\t" + node[0].getNumber() + "\t" + 847 node[1].getNumber() + "\t" + node[2].getNumber() + "\t" + 848 node[3].getNumber() + "\t" + node[4].getNumber() + "\t" + 849 node[5].getNumber() + "\t" + node[6].getNumber() + "\t" + 850 node[7].getNumber() + "\n" 851 ); 852 853 return out; 854 855 case RESULT_HEADER: 856 857 /* Print the header of the result file, for the block of Hexahedra (solid) elements. 858 * The element has eight gauss points where the results are calculated. They are placed on the natural positions 859 * which is where the gauss coordinates are most accurate. 860 */ 861 out = new String( 862 "GaussPoints \"Type" + type + 863 "\" ElemType Hexahedra \"MeshType" + type + "\"\n" 864 ); 865 out += ("Number Of Gauss Points: " + number_of_integration_points + 866 "\n"); 867 out += "Nodes Not Included\n"; // There are no gauss points in the nodes. 868 out += "Natural Coordinates: Given\n"; // They are on the optimum gauss coordinates instead, which GID will know by default when this switch is set to internal. 869 870 if (number_of_integration_points == 1) { 871 out += "0.0 0.0 0.0\n"; 872 } else if (number_of_integration_points == 8) { 873 out += "-0.57735 -0.57735 -0.57735\n"; 874 out += "0.57735 -0.57735 -0.57735\n"; 875 out += "0.57735 0.57735 -0.57735\n"; 876 out += "-0.57735 0.57735 -0.57735\n"; 877 out += "-0.57735 -0.57735 0.57735\n"; 878 out += "0.57735 -0.57735 0.57735\n"; 879 out += "0.57735 0.57735 0.57735\n"; 880 out += "-0.57735 0.57735 0.57735\n"; 881 } 882 883 out += "End GaussPoints\n"; 884 885 return out; 886 887 case RESULT_SUB_HEADER: 888 889 /* Print the subheader for the resultfile to initate each block of data from this element type 890 First parameter: Kind of results ( 1= scalar, 2=vector, 3=matrix, 4=2D plane deformation matrix, 5=Main stresser, 6=Euler angles (for local axes) 891 Second param: Location of the data (1= on the nodes, 2= in the gauss points); 892 Third param: Will there be a description on the results button in GID? (0= Make one up GID! 1= Yes, a description will be given) 893 Fourth param: Specify the name of the gauss point set that will be used "name" 894 */ 895 out = new String(" 3 2 0 \"Type" + type + "\"\n"); 896 897 return out; 898 899 case RESULT_STRESS_GLOBAL: 900 901 /* Print the Gauss stresses for this element and the requested gauss point. */ 902 if (gpn == 0) { 903 out = new String(number + " "); // Element number must start the first gauss point results 904 } else { 905 out = new String(""); // The rest of the lines will have no initial element number 906 } 907 908 /* 909 910 // The result is printed as [Sxx Syy Szz Sxy Syz Sxz] 911 out += (stress[gpn].get(0, 0) + "\t" + stress[gpn].get(1, 0) + 912 "\t" + stress[gpn].get(2, 0) + "\t" + stress[gpn].get(3, 0) + "\t" + 913 stress[gpn].get(4, 0) + "\t" + stress[gpn].get(5, 0) + "\n"); 914 */ 915 return out; 916 917 case RESULT_STRAIN_GLOBAL: 918 919 /* Print the Gauss stresses for this element and the requested gauss point. */ 920 if (gpn == 0) { 921 out = new String(number + " "); // Element number must start the first gauss point results 922 } else { 923 out = new String(""); // The rest of the lines will have no initial element number 924 } 925 926 927 /* 928 // The result is printed as [Sxx Syy Szz Sxy Syz Sxz] 929 out += (strain[gpn].get(0, 0) + "\t" + strain[gpn].get(1, 0) + 930 "\t" + strain[gpn].get(2, 0) + "\t" + strain[gpn].get(3, 0) + "\t" + 931 strain[gpn].get(4, 0) + "\t" + strain[gpn].get(5, 0) + "\n"); 932 */ 933 return out; 934 935 default: 936 return new String(""); 937 } 938 } 939 940 941 /** 942 * This method is used to create the lines needed in the result file. The 943 * method returns a string which is printed directly. However, due to the 944 * fact that the line may be different depending on what is requested to 945 * be printed and that the number of methods should be kept down, the 946 * first parameter here is a control parameter. This parameter describes 947 * what should be printed. The second parameter is a required input when 948 * gauss point results are to be printed. Creation date: (09/12/01 %T) 949 * 950 * @param ctrl The control number to say if a header of result file is to 951 * be printed. 952 * @param gpn The gauss point number. 953 * 954 * @return java.lang.String 955 */ print_Fembic(int ctrl, int gpn)956 public String print_Fembic(int ctrl, int gpn) { 957 String out; 958 959 switch (ctrl) { 960 961 case MESH: 962 963 /* Print the element number and connected nodes */ 964 out = new String( 965 number + "\t nodes = [" + node[0].getNumber() + "," 966 + node[1].getNumber() + "," 967 + node[2].getNumber() + "," 968 + node[3].getNumber() + "," 969 + node[4].getNumber() + "," 970 + node[5].getNumber() + "," 971 + node[6].getNumber() + "," 972 + node[7].getNumber() + "]\t" 973 + "material = " + material[0].getName() + "\t" 974 ); 975 976 977 if (NIP_is_set) 978 out += " nip = " + number_of_integration_points; 979 980 out += "\n"; 981 982 return out; 983 984 985 default: 986 return new String(""); 987 } 988 } 989 990 991 992 993 /** 994 * Insert the method's description here. Creation date: (27/09/01 %T) 995 */ setInitialConditions()996 public void setInitialConditions() 997 throws IllegalArgumentException 998 { 999 double t = 1.0 / Math.sqrt(3); 1000 int i; 1001 int j; 1002 1003 /* Make a local copy of the material object and keep it inside the element. 1004 * The reason for this is twofold. 1005 * 1. The material law is now tied to the element and will remember the element properties and history 1006 * 2. Since the material is cloned now, after Initialize, all the parameters are defined in the material object and are 1007 * now automatically transferred to the local copy. 1008 * 1009 * Note that one copy of the materiallaw is required for each intergration point in the element. 1010 */ 1011 try { 1012 for (i = 0; i < number_of_integration_points; i++) { 1013 material[i] = (Material) material[0].copy(); 1014 } 1015 } catch (CloneNotSupportedException e) { 1016 System.err.println("Object cannot clone"); 1017 } 1018 1019 // Now call any necessary initialisations in the law. 1020 for (i = 0; i < number_of_integration_points; i++) { 1021 material[i].setInitialConditions(); 1022 } 1023 1024 // 1025 // Define all important matrices 1026 // 1027 // For a single gauss point per element 1028 1029 /* 1030 if (number_of_integration_points == 1) { 1031 xsi[0] = 0; 1032 etha[0] = 0; 1033 phi[0] = 0; 1034 W[0] = 2.0; 1035 } else 1036 // For eight gauss points per element 1037 { 1038 xsi[0] = -t; 1039 xsi[1] = -t; 1040 xsi[2] = -t; 1041 xsi[3] = -t; 1042 xsi[4] = t; 1043 xsi[5] = t; 1044 xsi[6] = t; 1045 xsi[7] = t; 1046 1047 // 1048 etha[0] = -t; 1049 etha[1] = -t; 1050 etha[2] = t; 1051 etha[3] = t; 1052 etha[4] = -t; 1053 etha[5] = -t; 1054 etha[6] = t; 1055 etha[7] = t; 1056 1057 // 1058 phi[0] = t; 1059 phi[1] = -t; 1060 phi[2] = -t; 1061 phi[3] = t; 1062 phi[4] = t; 1063 phi[5] = -t; 1064 phi[6] = -t; 1065 phi[7] = t; 1066 1067 // 1068 W[0] = 1.0; 1069 W[1] = 1.0; 1070 W[2] = 1.0; 1071 W[3] = 1.0; 1072 W[4] = 1.0; 1073 W[5] = 1.0; 1074 W[6] = 1.0; 1075 W[7] = 1.0; 1076 } 1077 */ 1078 H.set(0, 0, 1.0); 1079 H.set(0, 1, 0); 1080 H.set(0, 2, 0); 1081 H.set(0, 3, 0); 1082 H.set(0, 4, 0); 1083 H.set(0, 5, 0); 1084 H.set(0, 6, 0); 1085 H.set(0, 7, 0); 1086 H.set(0, 8, 0); 1087 1088 // 1089 H.set(1, 0, 0); 1090 H.set(1, 1, 0); 1091 H.set(1, 2, 0); 1092 H.set(1, 3, 0); 1093 H.set(1, 4, 1.0); 1094 H.set(1, 5, 0); 1095 H.set(1, 6, 0); 1096 H.set(1, 7, 0); 1097 H.set(1, 8, 0); 1098 1099 // 1100 H.set(2, 0, 0); 1101 H.set(2, 1, 0); 1102 H.set(2, 2, 0); 1103 H.set(2, 3, 0); 1104 H.set(2, 4, 0); 1105 H.set(2, 5, 0); 1106 H.set(2, 6, 0); 1107 H.set(2, 7, 0); 1108 H.set(2, 8, 1.0); 1109 1110 // 1111 H.set(3, 0, 0); 1112 H.set(3, 1, 1.0); 1113 H.set(3, 2, 0); 1114 H.set(3, 3, 1.0); 1115 H.set(3, 4, 0); 1116 H.set(3, 5, 0); 1117 H.set(3, 6, 0); 1118 H.set(3, 7, 0); 1119 H.set(3, 8, 0); 1120 1121 // 1122 H.set(4, 0, 0); 1123 H.set(4, 1, 0); 1124 H.set(4, 2, 0); 1125 H.set(4, 3, 0); 1126 H.set(4, 4, 0); 1127 H.set(4, 5, 1.0); 1128 H.set(4, 6, 0); 1129 H.set(4, 7, 1.0); 1130 H.set(4, 8, 0); 1131 1132 // 1133 H.set(5, 0, 0); 1134 H.set(5, 1, 0); 1135 H.set(5, 2, 1.0); 1136 H.set(5, 3, 0); 1137 H.set(5, 4, 0); 1138 H.set(5, 5, 0); 1139 H.set(5, 6, 1.0); 1140 H.set(5, 7, 0); 1141 H.set(5, 8, 0); 1142 1143 // Set all elements in N to zero 1144 for (i = 0; i < 9; i++) { 1145 for (j = 0; j < 24; j++) { 1146 N.set(i, j, 0.0); 1147 } 1148 } 1149 1150 // Set all elements in P to zero 1151 for (i = 0; i < 9; i++) { 1152 for (j = 0; j < 9; j++) { 1153 P.set(i, j, 0.0); 1154 } 1155 } 1156 } 1157 1158 /** 1159 * This method calculates the local coordinate system for the element. The 1160 * method returns a handle to the system and this matrix is then stored in 1161 * the local_coordinate_system for later use. The matrix can be used in 1162 * the transformation of displacements and forces between the local and 1163 * global coordinates. However, in this element, the transformation is 1164 * made automatically in the matrix algebra of calculating the element 1165 * strains (they are derived in global directions) Creation date: 1166 * (25/12/01 %T) 1167 */ updateLocalCoordinateSystem()1168 public void updateLocalCoordinateSystem() { 1169 } 1170 1171 1172 } 1173 1174