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/Shell.java
23 import java.text.ParseException;
24
25 import java.util.*;
26
27
28 /**
29 * This is a Belytchko Tsai Lin element. It is based on material from various
30 * sources such as: Explicit Algorithms for the Nonlinear Dynamics of Shells -
31 * Computer Methods in Applied Mechanics and Engineering 42 (1984) 225-251
32 * Dyna Theoretical Manual Radioss Manual
33 *
34 * @author Jonas Forssell, Yuriy Mikhaylovskiy
35 *
36 * @see OtherClasses
37 */
38 public class Shell_BT_4 extends Element {
39 private static int DISABLED = 0;
40 private static int BASIC = 1;
41 private static int ADVANCED = 2;
42 private static int EDGE = 3;
43 private static int ADVANCED_EDGE = 4;
44 private double area;
45 private Jama.Matrix local_coordinate_system;
46 private Jama.Matrix[] strain;
47 private Jama.Matrix[] dstrain;
48 private Jama.Matrix strainrate;
49 private Jama.Matrix[] stress;
50 private Jama.Matrix force;
51 private Jama.Matrix hourglass_force;
52 private Jama.Matrix moment;
53 private Jama.Matrix hourglass_moment;
54 private Jama.Matrix B;
55 private Jama.Matrix p;
56 private Jama.Matrix v;
57 private Jama.Matrix a;
58 private Jama.Matrix trash;
59 private double dmx;
60 private double dmy;
61 private double dmxy;
62 private double kappax;
63 private double kappay;
64 private double kappaxy;
65 private double[] z;
66 private double[] weight_factor;
67 private int number_of_integration_points;
68 private int Contact;
69 private int printed_integration_point;
70 private double thickness;
71 private double factor;
72 private double friction;
73 private double original_thickness;
74 private Material[] material;
75 private double shear_factor;
76 private boolean hourglass_control;
77 private double membrane_hourglass_coeff;
78 private double out_of_plane_hourglass_coeff;
79 private double rotation_hourglass_coeff;
80 private Jama.Matrix hourglass_vector;
81 private Load load;
82 private boolean NIP_is_set;
83 private boolean PIP_is_set;
84 private boolean Nodes_are_set;
85 private boolean Material_is_set;
86 private boolean T_is_set;
87 private boolean Factor_is_set;
88 private boolean Friction_is_set;
89 private boolean Thinning_is_enabled;
90 private Contact_Triangle internal_contact_element_1;
91 private Contact_Triangle internal_contact_element_2;
92 private Contact_Triangle internal_contact_element_3;
93 private Contact_Triangle internal_contact_element_4;
94 private Contact_Line internal_contact_line_element_1;
95 private Contact_Line internal_contact_line_element_2;
96 private Contact_Line internal_contact_line_element_3;
97 private Contact_Line internal_contact_line_element_4;
98
99 // Attributes
100 // Associations
101 // Operations
102
103 /**
104 * Insert the method's description here. Creation date: (2001-08-10
105 * 19:45:45)
106 */
Shell_BT_4()107 public Shell_BT_4() {
108 super();
109 type = new String("SHELL_BT_4");
110
111 int i;
112 p = new Jama.Matrix(12, 1);
113 v = new Jama.Matrix(12, 1);
114 a = new Jama.Matrix(12, 1);
115 trash = new Jama.Matrix(12, 1);
116 B = new Jama.Matrix(2, 4);
117 node = new Node[4];
118 force = new Jama.Matrix(3, 1);
119 hourglass_force = new Jama.Matrix(3, 1);
120 moment = new Jama.Matrix(3, 1);
121 hourglass_moment = new Jama.Matrix(3, 1);
122 material = new Material[5];
123 strainrate = new Jama.Matrix(6, 1);
124 local_coordinate_system = new Jama.Matrix(3, 3);
125
126 //
127 strain = new Jama.Matrix[5];
128
129 for (i = 0; i < 5; i++) {
130 strain[i] = new Jama.Matrix(6, 1);
131 }
132
133 //
134 dstrain = new Jama.Matrix[5];
135
136 for (i = 0; i < 5; i++) {
137 dstrain[i] = new Jama.Matrix(6, 1);
138 }
139
140 //
141 stress = new Jama.Matrix[5];
142
143 for (i = 0; i < 5; i++) {
144 stress[i] = new Jama.Matrix(6, 1);
145 }
146
147 //
148 z = new double[5];
149 weight_factor = new double[5];
150
151 //
152 hourglass_vector = new Jama.Matrix(4, 1);
153
154 // Set the default values
155 // These will change if the user has defined them in the input file
156 shear_factor = 1.0;
157 hourglass_control = true;
158 membrane_hourglass_coeff = 0.1;
159 out_of_plane_hourglass_coeff = 0.1;
160 rotation_hourglass_coeff = 0.1;
161 Contact = BASIC;
162 Thinning_is_enabled = true;
163 number_of_integration_points = 5;
164 }
165
166 /**
167 * Insert the method's description here. Creation date: (25/12/01 %T)
168 */
assembleMassMatrix()169 public void assembleMassMatrix()
170 throws IllegalArgumentException
171 {
172 Jama.Matrix inertia;
173 inertia = new Jama.Matrix(3, 3);
174
175 double I;
176 double mass;
177
178 // Update local coordinate system
179 this.updateLocalCoordinateSystem();
180
181 // Determine the element Area
182 // and set the internal node coordinate (if applicable)
183 this.calculateLocalVariables();
184
185 // Distribute the mass contribution on the nodes
186 mass = material[0].getDensity() * area * thickness;
187 node[0].addMass(mass / 4);
188 node[1].addMass(mass / 4);
189 node[2].addMass(mass / 4);
190 node[3].addMass(mass / 4);
191
192 /* Now, calculate the inertia.
193 [Ixx Iyx Izx]
194 I = [Ixy Iyy Izy]
195 [Ixz Iyz Izz]
196 With the exact formula for intertia, the solution tends to diverge
197 for large rotation rates. A way to stabalize the rotation is to let
198 Izz = Iyy = Ixx , Ixy... = 0 and concider the rectangle as a square.
199 This gives the following:
200 */
201 I = (mass / 4) * ((area / 12) + ((thickness * thickness) / 12));
202
203 inertia.set(0, 0, I);
204 inertia.set(1, 1, I);
205 inertia.set(2, 2, I);
206 inertia.set(1, 0, 0);
207 inertia.set(2, 0, 0);
208 inertia.set(0, 1, 0);
209 inertia.set(0, 2, 0);
210 inertia.set(2, 1, 0);
211 inertia.set(1, 2, 0);
212
213 // Now, add to the nodes
214 node[0].addInertia(inertia);
215 node[1].addInertia(inertia);
216 node[2].addInertia(inertia);
217 node[3].addInertia(inertia);
218
219 //
220 if ((Contact == BASIC) || (Contact == EDGE)) {
221 internal_contact_element_1.assembleMassMatrix();
222 internal_contact_element_2.assembleMassMatrix();
223 }
224
225 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
226 internal_contact_element_1.assembleMassMatrix();
227 internal_contact_element_2.assembleMassMatrix();
228 internal_contact_element_3.assembleMassMatrix();
229 internal_contact_element_4.assembleMassMatrix();
230 }
231
232 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
233 if (internal_contact_line_element_1 != null) {
234 internal_contact_line_element_1.assembleMassMatrix();
235 }
236
237 if (internal_contact_line_element_2 != null) {
238 internal_contact_line_element_2.assembleMassMatrix();
239 }
240
241 if (internal_contact_line_element_3 != null) {
242 internal_contact_line_element_3.assembleMassMatrix();
243 }
244
245 if (internal_contact_line_element_4 != null) {
246 internal_contact_line_element_4.assembleMassMatrix();
247 }
248 }
249 }
250
251 /**
252 * Insert the method's description here. Creation date: (25/12/01 %T)
253 */
calculateContactForces()254 public void calculateContactForces() {
255 if ((Contact == BASIC) || (Contact == EDGE)) {
256 internal_contact_element_1.calculateContactForces();
257 internal_contact_element_2.calculateContactForces();
258 }
259
260 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
261 internal_contact_element_1.calculateContactForces();
262 internal_contact_element_2.calculateContactForces();
263 internal_contact_element_3.calculateContactForces();
264 internal_contact_element_4.calculateContactForces();
265
266 // Distribute the load from the internal node onto the others
267 node[0].addInternalForce(middle_node.getForce().times(0.25));
268 node[1].addInternalForce(middle_node.getForce().times(0.25));
269 node[2].addInternalForce(middle_node.getForce().times(0.25));
270 node[3].addInternalForce(middle_node.getForce().times(0.25));
271
272 // And reset it
273 middle_node.clearNodalForces();
274 }
275
276 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
277 if (internal_contact_line_element_1 != null) {
278 internal_contact_line_element_1.calculateContactForces();
279 }
280
281 if (internal_contact_line_element_2 != null) {
282 internal_contact_line_element_2.calculateContactForces();
283 }
284
285 if (internal_contact_line_element_3 != null) {
286 internal_contact_line_element_3.calculateContactForces();
287 }
288
289 if (internal_contact_line_element_4 != null) {
290 internal_contact_line_element_4.calculateContactForces();
291 }
292 }
293 }
294
295 /**
296 * Insert the method's description here. Creation date: (25/12/01 %T)
297 */
calculateExternalForces(double currtime)298 public void calculateExternalForces(double currtime) {
299 int i;
300
301 // Add the pressure force in local z-direction if there is a load object defined
302 if (load != null) {
303 force.set(0, 0, 0);
304 force.set(1, 0, 0);
305 force.set(2, 0, area * load.getPressure(currtime) * 0.25);
306
307 // Transform to global directions
308 force = local_coordinate_system.transpose().times(force);
309
310 // Distribute on nodes (note this is external forces so no sign reversal)
311 for (i = 0; i < 4; i++) {
312 node[i].addInternalForce(force);
313 }
314 }
315
316 if ((Contact == BASIC) || (Contact == EDGE)) {
317 internal_contact_element_1.calculateExternalForces(currtime);
318 internal_contact_element_2.calculateExternalForces(currtime);
319 }
320
321 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
322 internal_contact_element_1.calculateExternalForces(currtime);
323 internal_contact_element_2.calculateExternalForces(currtime);
324 internal_contact_element_3.calculateExternalForces(currtime);
325 internal_contact_element_4.calculateExternalForces(currtime);
326 }
327
328 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
329 if (internal_contact_line_element_1 != null) {
330 internal_contact_line_element_1.calculateExternalForces(
331 currtime
332 );
333 }
334
335 if (internal_contact_line_element_2 != null) {
336 internal_contact_line_element_2.calculateExternalForces(
337 currtime
338 );
339 }
340
341 if (internal_contact_line_element_3 != null) {
342 internal_contact_line_element_3.calculateExternalForces(
343 currtime
344 );
345 }
346
347 if (internal_contact_line_element_4 != null) {
348 internal_contact_line_element_4.calculateExternalForces(
349 currtime
350 );
351 }
352 }
353 }
354
355 /**
356 * This method is redesigned for the shell element. It calculates the local
357 * base vectors for the element. These vectors defines the local
358 * coordinate system for the element on which all the stress, strain and
359 * force calculations are based.
360 *
361 * The input are four points in space. They are used as follows:
362 * 1. The first point is the base point.
363 * 2. The local x-axis is defined from base point to the second point.
364 * 3. The local z-axis is defined as the cross product of the diagonals 1-3, 2-4
365 *
366 *
367 *
368 *
369 *
370 *
371 *
372 * 3. The local y-axis is defined from the base point to the third point 4. The cross
373 * vector between the x-axis and the y axis, defines the normal to the
374 * plane 5. A new adjusted x-axis is calculated. 6. Finally, the y-axis is
375 * calculated as the cross product of the x-axis and z-axis. Now, all the
376 * three axis are normal to each other and a correct coordinate system is
377 * defined. The method will return a matrix with these vectors: x-vec
378 * y-vec z-vec x y z Creation date:
379 * (26/08/01 %T) Jonas Forssell
380 */
calculateLocalBaseVectors( double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3, double x4, double y4, double z4, Jama.Matrix bvs )381 public void calculateLocalBaseVectors(
382 double x1, double y1, double z1, double x2, double y2, double z2,
383 double x3, double y3, double z3, double x4, double y4, double z4,
384 Jama.Matrix bvs
385 ) {
386 // 1&2. Define the local x-axis.
387 bvs.set(0, 0, x2 - x1);
388 bvs.set(0, 1, y2 - y1);
389 bvs.set(0, 2, z2 - z1);
390
391 // 3. Define the y-axis
392 bvs.set(1, 0, x4 - x1);
393 bvs.set(1, 1, y4 - y1);
394 bvs.set(1, 2, z4 - z1);
395
396 // 4. Calculate the z-axis (cross product of the diagonals)
397
398 bvs.set(2,0,(y3-y1)*(z4-z2)-(z3-z1)*(y4-y2));
399 bvs.set(2,1,(z3-z1)*(x4-x2)-(z4-z2)*(x3-x1));
400 bvs.set(2,2,(x3-x1)*(y4-y2)-(y3-y1)*(x4-x2));
401
402 // and normalize
403 //local_z_axis = local_z_axis.times(1.0/local_z_axis.length());
404 bvs.set(
405 1, 0,
406 Math.sqrt(
407 (bvs.get(2, 0) * bvs.get(2, 0)) +
408 (bvs.get(2, 1) * bvs.get(2, 1)) +
409 (bvs.get(2, 2) * bvs.get(2, 2))
410 )
411 );
412 bvs.set(2, 0, bvs.get(2, 0) / bvs.get(1, 0));
413 bvs.set(2, 1, bvs.get(2, 1) / bvs.get(1, 0));
414 bvs.set(2, 2, bvs.get(2, 2) / bvs.get(1, 0));
415
416 // 5. Adjust the x-axis (using the y-axis as a temporary variable)
417 // local_y_axis = local_x_axis.times(local_z_axis.transpose()).times(local_z_axis);
418 bvs.set(
419 1, 0,
420 (bvs.get(0, 0) * bvs.get(2, 0) * bvs.get(2, 0)) +
421 (bvs.get(0, 0) * bvs.get(2, 1) * bvs.get(2, 1)) +
422 (bvs.get(0, 0) * bvs.get(2, 2) * bvs.get(2, 2))
423 );
424 bvs.set(
425 1, 1,
426 (bvs.get(0, 1) * bvs.get(2, 0) * bvs.get(2, 0)) +
427 (bvs.get(0, 1) * bvs.get(2, 1) * bvs.get(2, 1)) +
428 (bvs.get(0, 1) * bvs.get(2, 2) * bvs.get(2, 2))
429 );
430 bvs.set(
431 1, 2,
432 (bvs.get(0, 2) * bvs.get(2, 0) * bvs.get(2, 0)) +
433 (bvs.get(0, 2) * bvs.get(2, 1) * bvs.get(2, 1)) +
434 (bvs.get(0, 2) * bvs.get(2, 2) * bvs.get(2, 2))
435 );
436
437 // local_x_axis = local_y_axis.times(1.0/local_y_axis.length());
438 bvs.set(
439 0, 2,
440 Math.sqrt(
441 (bvs.get(1, 0) * bvs.get(1, 0)) +
442 (bvs.get(1, 1) * bvs.get(1, 1)) +
443 (bvs.get(1, 2) * bvs.get(1, 2))
444 )
445 );
446 bvs.set(0, 0, bvs.get(1, 0) / bvs.get(0, 2));
447 bvs.set(0, 1, bvs.get(1, 1) / bvs.get(0, 2));
448 bvs.set(0, 2, bvs.get(1, 2) / bvs.get(0, 2));
449
450 // 6. Adjust the y-axis
451 // local_y_axis = local_z_axis.vectorProduct(local_x_axis);
452 bvs.set(
453 1, 0,
454 (bvs.get(2, 1) * bvs.get(0, 2)) - (bvs.get(2, 2) * bvs.get(0, 1))
455 );
456 bvs.set(
457 1, 1,
458 (bvs.get(2, 2) * bvs.get(0, 0)) - (bvs.get(2, 0) * bvs.get(0, 2))
459 );
460 bvs.set(
461 1, 2,
462 (bvs.get(2, 0) * bvs.get(0, 1)) - (bvs.get(2, 1) * bvs.get(0, 0))
463 );
464 }
465
466 /**
467 * This calculates and adds the element internal forces to the nodes. Note,
468 * by definition, the internal forces are subtracted from the external
469 * loads Thus, the "addition" is negative. Creation date: (2001-10-23
470 * 14.03.06)
471 *
472 */
calculateNodalForces(int integration_point, double timestep)473 public void calculateNodalForces(int integration_point, double timestep)
474 {
475 /* The nodal forces are computed as a function of the gauss integration of the
476 stresses. This means a summary must be made for all integration points.
477 This routine will be called once for every integration point so this routine
478 only has to calculate the force contribution of a specific integration point and
479 distribute it on the nodes. */
480 int i;
481 double ki;
482 double ko;
483 double kr;
484 Jama.Matrix transposed_coordinate_system;
485
486 ki = 0.0;
487 ko = 0.0;
488 kr = 0.0;
489 transposed_coordinate_system = local_coordinate_system.transpose();
490
491 // Check the area of the element
492 if (integration_point == 0) {
493 if (area <= 0) {
494 throw new IllegalArgumentException(
495 "Error in Shell_BT_4 element. Negative or zero area. "
496 );
497 }
498 }
499
500 // Loop through the nodes
501 for (i = 0; i < 4; i++) {
502 // Calculate forces fx, fy , fz
503 force.set(
504 0, 0,
505 0.5 * thickness * area * weight_factor[integration_point] * ((B.get(0,
506 i
507 ) * stress[integration_point].get(0, 0)) +
508 (B.get(1, i) * stress[integration_point].get(3, 0)))
509 );
510 force.set(
511 1, 0,
512 0.5 * thickness * area * weight_factor[integration_point] * ((B.get(1,
513 i
514 ) * stress[integration_point].get(1, 0)) +
515 (B.get(0, i) * stress[integration_point].get(3, 0)))
516 );
517 force.set(
518 2, 0,
519 0.5 * thickness * area * shear_factor * weight_factor[integration_point] * ((B.get(0,
520 i
521 ) * stress[integration_point].get(5, 0)) +
522 (B.get(1, i) * stress[integration_point].get(4, 0)))
523 );
524
525 /* Now, add the hourglass forces
526 These are calculated according to a viscous model where the forces are related to the
527 rate of displacement (node velocity). The formulation is based on work done by
528 Kosloff and Frasier.
529 */
530 if ((integration_point == 0) && (hourglass_control == true)) {
531 if (i == 0) {
532 ki = 0.25 * material[0].getDensity() * material[0].wavespeedTwoDimensional(
533 0.0, 0.0
534 ) * thickness * Math.sqrt(
535 (membrane_hourglass_coeff * area) / 2
536 );
537 ko = 0.25 * material[0].getDensity() * material[0].wavespeedTwoDimensional(
538 0.0, 0.0
539 ) * thickness * thickness * Math.sqrt(
540 out_of_plane_hourglass_coeff / 10
541 );
542 kr = 0.02 * material[0].getDensity() * material[0].wavespeedTwoDimensional(
543 0.0, 0.0
544 ) * thickness * thickness * area * Math.sqrt(
545 rotation_hourglass_coeff / 2
546 );
547 }
548
549 // In-plane mode (x and y directions)
550 hourglass_force.set(
551 0, 0,
552 ki * ((v.get(0, 0) - v.get(3, 0) + v.get(6, 0)) -
553 v.get(9, 0)) * hourglass_vector.get(i, 0)
554 );
555 hourglass_force.set(
556 1, 0,
557 ki * ((v.get(1, 0) - v.get(4, 0) + v.get(7, 0)) -
558 v.get(10, 0)) * hourglass_vector.get(i, 0)
559 );
560
561 // Out-of-plane mode (z directions)
562 hourglass_force.set(
563 2, 0,
564 ko * ((v.get(2, 0) - v.get(5, 0) + v.get(8, 0)) -
565 v.get(11, 0)) * hourglass_vector.get(i, 0)
566 );
567
568 // .. and transform to global coordinate system
569 hourglass_force = transposed_coordinate_system.times(
570 hourglass_force
571 );
572
573 // add to the node
574 node[i].addHourglassForce(hourglass_force.times(-1.0));
575 }
576
577 // .. and transform to global coordinate system
578 force = transposed_coordinate_system.times(force);
579
580 // Calculate moments mx, my, mz
581 // For those who are studying this in detail and is wondering why the formula is 0.5*thickness and not 0.25*thickness*thickness, remember that z[integration_point] already includes 0.5*thickness
582 moment.set(
583 0, 0,
584 0.5 * thickness * area * weight_factor[integration_point] * (((B.get(1,
585 i
586 ) * (-z[integration_point] * stress[integration_point].get(
587 1, 0
588 ))) +
589 (B.get(0, i) * (-z[integration_point] * stress[integration_point].get(
590 3, 0
591 )))) -
592 (0.25 * shear_factor * stress[integration_point].get(4, 0)))
593 );
594 moment.set(
595 1, 0,
596 0.5 * thickness * area * weight_factor[integration_point] * ((-B.get(0,
597 i
598 ) * (-z[integration_point] * stress[integration_point].get(
599 0, 0
600 ))) -
601 (B.get(1, i) * (-z[integration_point] * stress[integration_point].get(
602 3, 0
603 ))) +
604 (0.25 * shear_factor * stress[integration_point].get(5, 0)))
605 );
606 moment.set(2, 0, 0);
607
608 // Hourglass
609 if ((integration_point == 0) && (hourglass_control == true)) {
610 hourglass_moment.set(
611 0, 0,
612 kr * ((a.get(0, 0) - a.get(3, 0) + a.get(6, 0)) -
613 a.get(9, 0)) * hourglass_vector.get(i, 0)
614 );
615 hourglass_moment.set(
616 1, 0,
617 kr * ((a.get(1, 0) - a.get(4, 0) + a.get(7, 0)) -
618 a.get(10, 0)) * hourglass_vector.get(i, 0)
619 );
620 hourglass_moment.set(
621 2, 0,
622 kr * ((a.get(2, 0) - a.get(5, 0) + a.get(8, 0)) -
623 a.get(11, 0)) * hourglass_vector.get(i, 0)
624 );
625
626 // .. and transform to global coordinate system
627 hourglass_moment = transposed_coordinate_system.times(
628 hourglass_moment
629 );
630
631 // add to the node.
632 node[i].addHourglassMoment(hourglass_moment.times(-1.0));
633 }
634
635 // .. and transform to global coordinate system
636 moment = transposed_coordinate_system.times(moment);
637
638 /* Now, add the force and moment loads to the node. Remember that internal forces are
639 negative and external are positive. */
640 node[i].addInternalForce(force.times(-1.0));
641 node[i].addInternalMoment(moment.times(-1.0));
642 }
643 }
644
645 /**
646 * v is an array of node velocities expressed in the local coordinate
647 * system v = [u1 v1 w1 u2 v2 w2 u3 v3 w3 u4 v4 w4] p is an array of node
648 * coordinates expessed in the local coordinate system p = [x1 y1 z1 x2 y2
649 * z2 x3 y3 z3 x4 y4 z4] a is an array of node rotation velocity expressed
650 * in the local coordinate system a = [ax1 ay1 az1 ax2 ay2 az2 ax3 ay3 az3
651 * ax4 ay4 az4] B is a matrix of the shape function derivatives. Creation
652 * date: (25/12/01 %T) Jonas Forssell
653 *
654 * @param tstep double
655 */
calculateStrain(double tstep, int integration_point)656 public void calculateStrain(double tstep, int integration_point) {
657 // Do this block only once for each timestep. Not for each integration point
658 if (integration_point == 0) {
659 this.updateLocalCoordinateSystem();
660 this.calculateLocalVariables();
661
662 // Correct the area if element is set to high accuracy
663 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
664 area = internal_contact_element_1.getArea() +
665 internal_contact_element_2.getArea() +
666 internal_contact_element_3.getArea() +
667 internal_contact_element_4.getArea();
668 }
669
670 /* Now that we have the input in local coordinate system, it
671 is time to calculate the strain rates.
672 Start by creating the B matrix (the derivate of the shape functions)
673 which will be used frequently in the rest of the calculations. */
674 B.set(0, 0, p.get(4, 0) - p.get(10, 0));
675 B.set(0, 1, p.get(7, 0) - p.get(1, 0));
676 B.set(0, 2, p.get(10, 0) - p.get(4, 0));
677 B.set(0, 3, p.get(1, 0) - p.get(7, 0));
678 B.set(1, 0, p.get(9, 0) - p.get(3, 0));
679 B.set(1, 1, p.get(0, 0) - p.get(6, 0));
680 B.set(1, 2, p.get(3, 0) - p.get(9, 0));
681 B.set(1, 3, p.get(6, 0) - p.get(0, 0));
682 B.timesEquals(1.0 / (2.0 * area));
683
684 /* Now, the strain rates can be determined. They are calculated in
685 the middle of the shell and then transformed out through the thickness.
686
687 Start by calculating the middle strain rates: */
688 dmx = (B.get(0, 0) * (v.get(0, 0) - v.get(6, 0))) +
689 (B.get(0, 1) * (v.get(3, 0) - v.get(9, 0)));
690 dmy = (B.get(1, 0) * (v.get(1, 0) - v.get(7, 0))) +
691 (B.get(1, 1) * (v.get(4, 0) - v.get(10, 0)));
692 dmxy = (B.get(1, 0) * (v.get(0, 0) - v.get(6, 0))) +
693 (B.get(1, 1) * (v.get(3, 0) - v.get(9, 0))) +
694 (B.get(0, 0) * (v.get(1, 0) - v.get(7, 0))) +
695 (B.get(0, 1) * (v.get(4, 0) - v.get(10, 0)));
696
697 //
698 kappax = -((B.get(0, 0) * (a.get(1, 0) - a.get(7, 0))) +
699 (B.get(0, 1) * (a.get(4, 0) - a.get(10, 0))));
700 kappay = (B.get(1, 0) * (a.get(0, 0) - a.get(6, 0))) +
701 (B.get(1, 1) * (a.get(3, 0) - a.get(9, 0)));
702 kappaxy = (-B.get(1, 0) * (a.get(1, 0) - a.get(7, 0))) -
703 (B.get(1, 1) * (a.get(4, 0) - a.get(10, 0))) +
704 (B.get(0, 0) * (a.get(0, 0) - a.get(6, 0))) +
705 (B.get(0, 1) * (a.get(3, 0) - a.get(9, 0)));
706 }
707
708 // Now, compute the strainrate for this integration point.
709 strainrate.set(0, 0, dmx - (z[integration_point] * kappax));
710 strainrate.set(1, 0, dmy - (z[integration_point] * kappay));
711 strainrate.set(3, 0, dmxy - (z[integration_point] * kappaxy));
712
713 // The bending strain rates
714 // dyz
715 strainrate.set(
716 4, 0,
717 ((B.get(1, 0) * (v.get(2, 0) - v.get(8, 0))) +
718 (B.get(1, 1) * (v.get(5, 0) - v.get(11, 0)))) -
719 (0.25 * (a.get(0, 0) + a.get(3, 0) + a.get(6, 0) + a.get(9, 0)))
720 );
721
722 // dxz
723 strainrate.set(
724 5, 0,
725 (B.get(0, 0) * (v.get(2, 0) - v.get(8, 0))) +
726 (B.get(0, 1) * (v.get(5, 0) - v.get(11, 0))) +
727 (0.25 * (a.get(1, 0) + a.get(4, 0) + a.get(7, 0) + a.get(10, 0)))
728 );
729
730 // Knowing the strainrate, the strain increase for this integration point is easily computed.
731 dstrain[integration_point] = strainrate.times(tstep);
732 }
733
734 /**
735 * Insert the method's description here. Creation date: (25/12/01 %T)
736 */
calculateStress(int integration_point, double timestep)737 public void calculateStress(int integration_point, double timestep) {
738 /* This calculates the stresses using the strain and strain increase.
739 The result is put back into the stress matrix. The element contains one local
740 material object for each integration point since some materials has a memory which must
741 be maintaned. */
742 material[integration_point].calculateStressTwoDimensionalPlaneStress(
743 strain[integration_point], dstrain[integration_point],
744 stress[integration_point], timestep
745 );
746
747 // Update the thickness for thickness reduction using the mid-most integration point.
748 if (
749 Thinning_is_enabled &&
750 (integration_point == (number_of_integration_points / 2))
751 ) {
752 thickness = (1 +
753 strain[number_of_integration_points / 2].get(2, 0)) * original_thickness;
754 }
755 }
756
757 /**
758 * Insert the method's description here. Creation date: (25/12/01 %T)
759 *
760 * @param current_timestep double
761 *
762 * @return double
763 */
checkTimestep(double current_timestep)764 public double checkTimestep(double current_timestep) {
765 double critical_length;
766 double timestep;
767
768 // Find the longest edge on the element
769 critical_length = Math.sqrt(
770 ((node[1].getX_pos() - node[0].getX_pos()) * (node[1].getX_pos() -
771 node[0].getX_pos())) +
772 ((node[1].getY_pos() - node[0].getY_pos()) * (node[1].getY_pos() -
773 node[0].getY_pos())) +
774 ((node[1].getZ_pos() - node[0].getZ_pos()) * (node[1].getZ_pos() -
775 node[0].getZ_pos()))
776 );
777 critical_length = Math.max(
778 critical_length,
779 Math.sqrt(
780 ((node[2].getX_pos() - node[1].getX_pos()) * (node[2].getX_pos() -
781 node[1].getX_pos())) +
782 ((node[2].getY_pos() - node[1].getY_pos()) * (node[2].getY_pos() -
783 node[1].getY_pos())) +
784 ((node[2].getZ_pos() - node[1].getZ_pos()) * (node[2].getZ_pos() -
785 node[1].getZ_pos()))
786 )
787 );
788 critical_length = Math.max(
789 critical_length,
790 Math.sqrt(
791 ((node[3].getX_pos() - node[2].getX_pos()) * (node[3].getX_pos() -
792 node[2].getX_pos())) +
793 ((node[3].getY_pos() - node[2].getY_pos()) * (node[3].getY_pos() -
794 node[2].getY_pos())) +
795 ((node[3].getZ_pos() - node[2].getZ_pos()) * (node[3].getZ_pos() -
796 node[2].getZ_pos()))
797 )
798 );
799 critical_length = Math.max(
800 critical_length,
801 Math.sqrt(
802 ((node[0].getX_pos() - node[3].getX_pos()) * (node[0].getX_pos() -
803 node[3].getX_pos())) +
804 ((node[0].getY_pos() - node[3].getY_pos()) * (node[0].getY_pos() -
805 node[3].getY_pos())) +
806 ((node[0].getZ_pos() - node[3].getZ_pos()) * (node[0].getZ_pos() -
807 node[3].getZ_pos()))
808 )
809 );
810 critical_length = area / critical_length;
811
812 // Determine the timestep.
813 timestep = critical_length / material[0].wavespeedTwoDimensional(
814 0.0, 0.0
815 );
816
817 // Now, return the smallest of the suggested and the calculated timestep.
818 // The smallest element will decide the timestep for the whole simulation.
819 return Math.min(0.96 * timestep, current_timestep);
820 }
821
822 /**
823 * This element only has one integration point (gauss point) Creation date:
824 * (26/12/01 %T)
825 *
826 * @return int
827 */
getNumberOfIntegrationPoints()828 public int getNumberOfIntegrationPoints() {
829 return number_of_integration_points;
830 }
831
832 /**
833 * Insert the method's description here. Creation date: (25/12/01 %T)
834 *
835 * @param arg1 java.lang.String
836 * @param arg2 java.lang.String
837 * @param arg3 java.lang.String
838 * @param lineno int
839 * @param nodelist java.util.Vector
840 * @param materiallist java.util.Vector
841 */
parse_Fembic( Token[] param, int lineno, RplVector nodelist, RplVector materiallist, RplVector loadlist, Hashtable nodetable )842 public void parse_Fembic(
843 Token[] param, int lineno, RplVector nodelist, RplVector materiallist,
844 RplVector loadlist, Hashtable nodetable
845 )
846 throws java.text.ParseException
847 {
848 int j;
849 int i = 0;
850 int index;
851 Load temp_load;
852 Token[] contact_input;
853
854 while (i < param.length) {
855 // The nodes of the element are defined
856 if (
857 param[i].getw().toUpperCase().equals("NODES") &&
858 param[i + 1].getw().toUpperCase().equals("=")
859 ) {
860 // Assume now that the nodes are delivered in param3, with the format
861 // [nodenr,nodenr,nodenr,nodenr]
862 if (
863 ! param[i + 2].getw().toUpperCase().startsWith("[") ||
864 ! param[i + 2].getw().toUpperCase().endsWith("]")
865 ) {
866 throw new java.text.ParseException(
867 "Error, node number definition should be [nodenr1,nodenr2,nodenr3,nodenr4]",
868 lineno
869 );
870 }
871
872 // Ok, now find the numbers
873 try {
874 for (j = 0; j < 4; j++) {
875 node[j] = super.findNode(
876 super.getNodeNumber(j + 1, param[i + 2].getw()),
877 nodetable
878 );
879 }
880
881 i += 3;
882 Nodes_are_set = true;
883 } catch (IllegalArgumentException e) {
884 throw new ParseException(e.getMessage(), lineno);
885 }
886 } else
887 // The material of the element is defined
888 if (
889 param[i].getw().toUpperCase().equals("MATERIAL") &&
890 param[i + 1].getw().toUpperCase().equals("=")
891 ) {
892 // Assume now that the material name is delivered in param3
893 try {
894 // We want the handle to the material object.
895 material[0] = super.findMaterial(
896 param[i + 2].getw().toUpperCase(), materiallist
897 );
898 i += 3;
899 Material_is_set = true;
900 } catch (IllegalArgumentException e) {
901 throw new java.text.ParseException(
902 "Error in Shell_BT_4 element\n" + e.getMessage(), lineno
903 );
904 }
905 } else
906 // The number of integration points of the element is defined
907 if (
908 param[i].getw().toUpperCase().equals("NIP") &&
909 param[i + 1].getw().toUpperCase().equals("=")
910 ) {
911 // The value of the cross section area is in param3. Set this in the element
912 number_of_integration_points = (int) param[i + 2].getn();
913 i += 3;
914 NIP_is_set = true;
915 } else
916 // The number of the integration point to print results from
917 if (
918 param[i].getw().toUpperCase().equals("PIP") &&
919 param[i + 1].getw().toUpperCase().equals("=")
920 ) {
921 // The value of the cross section area is in param3. Set this in the element
922 printed_integration_point = (int) param[i + 2].getn() - 1;
923 i += 3;
924 PIP_is_set = true;
925 } else
926 // The thickness of the shell is defined
927 if (
928 param[i].getw().toUpperCase().equals("T") &&
929 param[i + 1].getw().toUpperCase().equals("=")
930 ) {
931 // The value of the element thickness is in param3. Set this in the element
932 thickness = param[i + 2].getn();
933 original_thickness = thickness;
934 i += 3;
935 T_is_set = true;
936 } else
937 // The friction of the shell is defined
938 if (
939 param[i].getw().toUpperCase().equals("FRICTION") &&
940 param[i + 1].getw().toUpperCase().equals("=")
941 ) {
942 // The value of the element thickness is in param3. Set this in the element
943 friction = param[i + 2].getn();
944 i += 3;
945 Friction_is_set = true;
946 } else
947 // The shear factor of the shell is defined
948 if (
949 param[i].getw().toUpperCase().equals("SHEAR_FACTOR") &&
950 param[i + 1].getw().toUpperCase().equals("=")
951 ) {
952 shear_factor = param[i + 2].getn();
953 i += 3;
954 } else
955 // Hourglass control of the shell element is defined
956 // This can be either ON or OFF
957 if (
958 param[i].getw().toUpperCase().equals("HOURGLASS") &&
959 param[i + 1].getw().toUpperCase().equals("=")
960 ) {
961 // Default is horglass control is ON. Anything other than OFF means default.
962 if (param[i + 2].getw().toUpperCase().equals("OFF")) {
963 hourglass_control = false;
964 }
965
966 i += 3;
967 } else
968 // The shell membrane hourglass coefficient
969 if (
970 param[i].getw().toUpperCase().equals("MHC") &&
971 param[i + 1].getw().toUpperCase().equals("=")
972 ) {
973 membrane_hourglass_coeff = param[i + 2].getn();
974 i += 3;
975 } else
976 // The shell out of plane hourglass coefficient
977 if (
978 param[i].getw().toUpperCase().equals("OOPHC") &&
979 param[i + 1].getw().toUpperCase().equals("=")
980 ) {
981 out_of_plane_hourglass_coeff = param[i + 2].getn();
982 i += 3;
983 } else
984 // The shell rotation hourglass coefficient
985 if (
986 param[i].getw().toUpperCase().equals("RHC") &&
987 param[i + 1].getw().toUpperCase().equals("=")
988 ) {
989 rotation_hourglass_coeff = param[i + 2].getn();
990 i += 3;
991 } else if (
992 param[i].getw().toUpperCase().equals("LOAD") &&
993 param[i + 2].is_a_word()
994 ) {
995 for (index = 0; index < loadlist.size(); index++) {
996 temp_load = (Load) loadlist.elementAt(index);
997
998 if (
999 temp_load.name.equals(
1000 param[i + 2].getw().toUpperCase()
1001 )
1002 ) {
1003 load = temp_load;
1004
1005 break;
1006 }
1007 }
1008
1009 if (index == loadlist.size()) {
1010 throw new java.text.ParseException(
1011 "Load name specified does not exist", lineno
1012 );
1013 }
1014
1015 i += 3;
1016 } else
1017 // The factor for the internal contact element is defined
1018 if (
1019 param[i].getw().toUpperCase().equals("FACTOR") &&
1020 param[i + 1].getw().toUpperCase().equals("=")
1021 ) {
1022 // The value of the element thickness is in param3. Set this in the element
1023 factor = param[i + 2].getn();
1024 i += 3;
1025 Factor_is_set = true;
1026 } else
1027 // The contact element of the shell is specified
1028 if (
1029 param[i].getw().toUpperCase().equals("CONTACT") &&
1030 param[i + 1].getw().toUpperCase().equals("=")
1031 ) {
1032 // The value of the element thickness is in param3. Set this in the element
1033 if (param[i + 2].getw().toUpperCase().equals("OFF")) {
1034 Contact = DISABLED;
1035 } else if (param[i + 2].getw().toUpperCase().equals("ADVANCED")) {
1036 Contact = ADVANCED;
1037 } else if (param[i + 2].getw().toUpperCase().equals("EDGE")) {
1038 Contact = EDGE;
1039 } else if (
1040 param[i + 2].getw().toUpperCase().equals("ADVANCED_EDGE")
1041 ) {
1042 Contact = ADVANCED_EDGE;
1043 } else {
1044 throw new ParseException(
1045 "Unrecognized contact parameter", lineno
1046 );
1047 }
1048
1049 i += 3;
1050 } else
1051 // The thickness reduction feature of the shell is specified
1052 if (
1053 param[i].getw().toUpperCase().equals("THINNING") &&
1054 param[i + 1].getw().toUpperCase().equals("=")
1055 ) {
1056 // The value of the parameter is in param3.
1057 if (param[i + 2].getw().toUpperCase().equals("OFF")) {
1058 Thinning_is_enabled = false;
1059 } else {
1060 throw new ParseException(
1061 "Unrecognized thinning parameter", lineno
1062 );
1063 }
1064
1065 i += 3;
1066 } else {
1067 // Neither material or nodes are defined. Then the parameter is wrong.
1068 throw new java.text.ParseException(
1069 "Unknown Shell element parameter ", lineno
1070 );
1071 }
1072 }
1073
1074 // **** Time to parse and initialize the contact elements *****
1075 // Make the input for the contact element 1
1076 i = 6;
1077
1078 if (Factor_is_set) {
1079 i += 3;
1080 }
1081
1082 if (Friction_is_set) {
1083 i += 3;
1084 }
1085
1086 contact_input = new Token[i];
1087
1088 // Make a local instance of a contact element to handle the contacts
1089 i = 6;
1090
1091 if ((Contact == BASIC) || (Contact == EDGE)) {
1092 internal_contact_element_1 = new Contact_Triangle();
1093 internal_contact_element_2 = new Contact_Triangle();
1094
1095 internal_contact_element_1.setNumber(-1);
1096 internal_contact_element_2.setNumber(-2);
1097
1098 contact_input[0] = new Token(new String("nodes"));
1099 contact_input[1] = new Token(new String("="));
1100 contact_input[2] = new Token(
1101 new String(
1102 "[" + node[0].getNumber() + "," + node[1].getNumber() +
1103 "," + node[2].getNumber() + "]"
1104 )
1105 );
1106 contact_input[3] = new Token(new String("t"));
1107 contact_input[4] = new Token(new String("="));
1108 contact_input[5] = new Token(thickness);
1109
1110 if (Factor_is_set) {
1111 contact_input[i] = new Token(new String("factor"));
1112 contact_input[i + 1] = new Token(new String("="));
1113 contact_input[i + 2] = new Token(factor);
1114 i += 3;
1115 }
1116
1117 if (Friction_is_set) {
1118 contact_input[i] = new Token(new String("friction"));
1119 contact_input[i + 1] = new Token(new String("="));
1120 contact_input[i + 2] = new Token(friction);
1121 i += 3;
1122 }
1123
1124 // Now, set the input for the contact element
1125 internal_contact_element_1.parse_Fembic(
1126 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1127 );
1128
1129 // Make the input for the contact element 2
1130 contact_input[2] = new Token(
1131 new String(
1132 "[" + node[2].getNumber() + "," + node[3].getNumber() +
1133 "," + node[0].getNumber() + "]"
1134 )
1135 );
1136
1137 // Now, set the input for the contact element also
1138 internal_contact_element_2.parse_Fembic(
1139 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1140 );
1141 }
1142
1143 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
1144 internal_contact_element_1 = new Contact_Triangle();
1145 internal_contact_element_2 = new Contact_Triangle();
1146 internal_contact_element_3 = new Contact_Triangle();
1147 internal_contact_element_4 = new Contact_Triangle();
1148
1149 internal_contact_element_1.setNumber(-1);
1150 internal_contact_element_2.setNumber(-2);
1151 internal_contact_element_3.setNumber(-3);
1152 internal_contact_element_4.setNumber(-4);
1153
1154 // Set up a middle node to connect all contact triangles to
1155 // Note that it is an internal node which makes it special compared to a
1156 // standard node.
1157 // Use a negative number which is same as element number to make it unique
1158
1159 middle_node = new InternalNode();
1160 middle_node.setNumber(-number);
1161 middle_node.registerMasterElement(this);
1162
1163 // Add the internal node to the nodelist so it becomes part of the main loop
1164 nodelist.add(middle_node);
1165 nodetable.put(new Integer(-number),middle_node);
1166
1167 contact_input[0] = new Token(new String("nodes"));
1168 contact_input[1] = new Token(new String("="));
1169 contact_input[2] = new Token(
1170 new String(
1171 "[" + node[0].getNumber() + "," + node[1].getNumber() +
1172 "," + middle_node.getNumber() + "]"
1173 )
1174 );
1175 contact_input[3] = new Token(new String("t"));
1176 contact_input[4] = new Token(new String("="));
1177 contact_input[5] = new Token(thickness);
1178
1179 if (Factor_is_set) {
1180 contact_input[i] = new Token(new String("factor"));
1181 contact_input[i + 1] = new Token(new String("="));
1182 contact_input[i + 2] = new Token(factor);
1183 i += 3;
1184 }
1185
1186 if (Friction_is_set) {
1187 contact_input[i] = new Token(new String("friction"));
1188 contact_input[i + 1] = new Token(new String("="));
1189 contact_input[i + 2] = new Token(friction);
1190 i += 3;
1191 }
1192
1193 // Now, set the input for the contact element
1194 internal_contact_element_1.parse_Fembic(
1195 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1196 );
1197
1198 // Make the input for the contact element 2
1199 contact_input[2] = new Token(
1200 new String(
1201 "[" + node[1].getNumber() + "," + node[2].getNumber() +
1202 "," + middle_node.getNumber() + "]"
1203 )
1204 );
1205 internal_contact_element_2.parse_Fembic(
1206 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1207 );
1208
1209 // Make the input for the contact element 3
1210 contact_input[2] = new Token(
1211 new String(
1212 "[" + node[2].getNumber() + "," + node[3].getNumber() +
1213 "," + middle_node.getNumber() + "]"
1214 )
1215 );
1216 internal_contact_element_3.parse_Fembic(
1217 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1218 );
1219
1220 // Make the input for the contact element 4
1221 contact_input[2] = new Token(
1222 new String(
1223 "[" + node[3].getNumber() + "," + node[0].getNumber() +
1224 "," + middle_node.getNumber() + "]"
1225 )
1226 );
1227 internal_contact_element_4.parse_Fembic(
1228 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1229 );
1230
1231 }
1232
1233 /* Make local instances of the edges contact elements to handle the contacts.
1234 * Note that some of the contact_inputs are assumed set in the previous block */
1235 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
1236 if (! node[0].hasContact_LineElementConnectedTo(node[1])) {
1237 internal_contact_line_element_1 = new Contact_Line();
1238 internal_contact_line_element_1.setNumber(-1);
1239
1240 contact_input[2] = new Token(
1241 new String(
1242 "[" + node[0].getNumber() + "," +
1243 node[1].getNumber() + "]"
1244 )
1245 );
1246 contact_input[3] = new Token(new String("D"));
1247 contact_input[4] = new Token(new String("="));
1248 contact_input[5] = new Token(thickness);
1249
1250 // Set the input
1251 internal_contact_line_element_1.parse_Fembic(
1252 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1253 );
1254 }
1255
1256 if (! node[1].hasContact_LineElementConnectedTo(node[2])) {
1257 internal_contact_line_element_2 = new Contact_Line();
1258 internal_contact_line_element_2.setNumber(-2);
1259
1260 contact_input[2] = new Token(
1261 new String(
1262 "[" + node[1].getNumber() + "," +
1263 node[2].getNumber() + "]"
1264 )
1265 );
1266 contact_input[3] = new Token(new String("D"));
1267 contact_input[4] = new Token(new String("="));
1268 contact_input[5] = new Token(thickness);
1269
1270 // Set the input
1271 internal_contact_line_element_2.parse_Fembic(
1272 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1273 );
1274 }
1275
1276 if (! node[2].hasContact_LineElementConnectedTo(node[3])) {
1277 internal_contact_line_element_3 = new Contact_Line();
1278 internal_contact_line_element_3.setNumber(-3);
1279
1280 contact_input[2] = new Token(
1281 new String(
1282 "[" + node[2].getNumber() + "," +
1283 node[3].getNumber() + "]"
1284 )
1285 );
1286 contact_input[3] = new Token(new String("D"));
1287 contact_input[4] = new Token(new String("="));
1288 contact_input[5] = new Token(thickness);
1289
1290 // Set the input
1291 internal_contact_line_element_3.parse_Fembic(
1292 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1293 );
1294 }
1295
1296 if (! node[3].hasContact_LineElementConnectedTo(node[0])) {
1297 internal_contact_line_element_4 = new Contact_Line();
1298 internal_contact_line_element_4.setNumber(-4);
1299
1300 contact_input[2] = new Token(
1301 new String(
1302 "[" + node[3].getNumber() + "," +
1303 node[0].getNumber() + "]"
1304 )
1305 );
1306 contact_input[3] = new Token(new String("D"));
1307 contact_input[4] = new Token(new String("="));
1308 contact_input[5] = new Token(thickness);
1309
1310 // Set the input
1311 internal_contact_line_element_4.parse_Fembic(
1312 contact_input, lineno, nodelist, materiallist, loadlist, nodetable
1313 );
1314 }
1315 }
1316 }
1317
1318
1319 /**
1320 * Checks that all mandatory variables have been set
1321 */
checkIndata()1322 public void checkIndata()
1323 throws IllegalArgumentException
1324 {
1325 // Check that all the required parameters have been parsed
1326 if (! Nodes_are_set) {
1327 throw new IllegalArgumentException(
1328 "No nodes defined for Shell_BT_4 element nr" + number
1329 );
1330 }
1331
1332 if (! Material_is_set) {
1333 throw new IllegalArgumentException(
1334 "No Material defined for Shell_BT_4 element nr" + number
1335 );
1336 }
1337
1338 if (! T_is_set) {
1339 throw new IllegalArgumentException(
1340 "No Thickness (T) defined for Shell_BT_4 element nr" + number
1341 );
1342 }
1343
1344 if (PIP_is_set) {
1345 if (printed_integration_point >= number_of_integration_points) {
1346 throw new IllegalArgumentException(
1347 "Printed integration point larger than available points"
1348 );
1349 }
1350
1351 if (printed_integration_point < 0) {
1352 throw new IllegalArgumentException("Printed integration point less than 1");
1353 }
1354 }
1355
1356 //
1357 if ((Contact == BASIC) || (Contact == EDGE)) {
1358 internal_contact_element_1.checkIndata();
1359 internal_contact_element_2.checkIndata();
1360 }
1361
1362 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
1363 internal_contact_element_1.checkIndata();
1364 internal_contact_element_2.checkIndata();
1365 internal_contact_element_3.checkIndata();
1366 internal_contact_element_4.checkIndata();
1367 }
1368
1369 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
1370 if (internal_contact_line_element_1 != null) {
1371 internal_contact_line_element_1.checkIndata();
1372 }
1373
1374 if (internal_contact_line_element_2 != null) {
1375 internal_contact_line_element_2.checkIndata();
1376 }
1377
1378 if (internal_contact_line_element_3 != null) {
1379 internal_contact_line_element_3.checkIndata();
1380 }
1381
1382 if (internal_contact_line_element_4 != null) {
1383 internal_contact_line_element_4.checkIndata();
1384 }
1385 }
1386 }
1387
1388 /**
1389 *
1390 * The method checks destroyed element or not.
1391 * If the element is destroyed variable failed = true else failed = false.
1392 *
1393 */
checkIfFailed()1394 public void checkIfFailed() {
1395 if (
1396 ! material[0].failureStrainIsSet() &&
1397 ! material[0].failureStressIsSet()
1398 ) {
1399 failed = false;
1400 return;
1401 }
1402
1403 if (material[0].failureStressIsSet()) {
1404 //0.70710678118654752440084436210485 = Math.sqrt(2) / 2.0
1405 double s = 0.7071f * Math.sqrt(
1406 Math.pow(
1407 stress[printed_integration_point].get(0, 0) -
1408 stress[printed_integration_point].get(1, 0), 2
1409 ) +
1410 Math.pow(
1411 stress[printed_integration_point].get(1, 0) -
1412 stress[printed_integration_point].get(2, 0), 2
1413 ) +
1414 Math.pow(
1415 stress[printed_integration_point].get(2, 0) -
1416 stress[printed_integration_point].get(0, 0), 2
1417 ) +
1418 (6 * (Math.pow(
1419 stress[printed_integration_point].get(3, 0), 2
1420 ) +
1421 Math.pow(stress[printed_integration_point].get(4, 0), 2) +
1422 Math.pow(stress[printed_integration_point].get(5, 0), 2)))
1423 );
1424
1425 if (s > material[0].getFailureStress()) {
1426 failed = true;
1427 return;
1428 }
1429 }
1430
1431 if (material[0].failureStrainIsSet()) {
1432 //0.47140452079103168293389624140323 = Math.sqrt(2) / 3.0
1433 double e = 0.4714f * Math.sqrt(
1434 Math.pow(
1435 strain[printed_integration_point].get(0, 0) -
1436 strain[printed_integration_point].get(1, 0), 2
1437 ) +
1438 Math.pow(
1439 strain[printed_integration_point].get(1, 0) -
1440 strain[printed_integration_point].get(2, 0), 2
1441 ) +
1442 Math.pow(
1443 strain[printed_integration_point].get(2, 0) -
1444 strain[printed_integration_point].get(0, 0), 2
1445 ) +
1446 (1.5 * (Math.pow(
1447 strain[printed_integration_point].get(3, 0), 2
1448 ) +
1449 Math.pow(strain[printed_integration_point].get(4, 0), 2) +
1450 Math.pow(strain[printed_integration_point].get(5, 0), 2)))
1451 );
1452
1453 if (e > material[0].getFailureStrain()) {
1454 failed = true;
1455 return;
1456 }
1457 }
1458
1459 failed = false;
1460 }
1461
1462 /**
1463 * This method is used to create the lines needed in the result file. The
1464 * method returns a string which is printed directly. However, due to the
1465 * fact that the line may be different depending on what is requested to
1466 * be printed and that the number of methods should be kept down, the
1467 * first parameter here is a control parameter. This parameter describes
1468 * what should be printed. The second parameter is a required input when
1469 * gauss point results are to be printed. Creation date: (09/12/01 %T)
1470 *
1471 * @param ctrl int
1472 * @param gpn int
1473 *
1474 * @return java.lang.String
1475 */
print_Gid(int ctrl, int gpn)1476 public String print_Gid(int ctrl, int gpn) {
1477 String out;
1478 int i;
1479 switch (ctrl) {
1480 case MESH_HEADER:
1481
1482 /* Print the header for the mesh.
1483 * In this case, the type of element is Quadilateral and it uses 4 nodes
1484 */
1485 out = new String(
1486 "MESH \"MeshType" + type +
1487 "\" Dimension 3 ElemType Quadrilateral Nnode 4\n"
1488 );
1489
1490 return out;
1491
1492 case MESH:
1493
1494 /* Print the element number and connected nodes */
1495 out = new String(
1496 number + "\t" + node[0].getNumber() + "\t" +
1497 node[1].getNumber() + "\t" + node[2].getNumber() + "\t" +
1498 node[3].getNumber() + "\n"
1499 );
1500
1501 return out;
1502
1503 case RESULT_HEADER:
1504
1505 /* Print the header of the result file, for the block of Quadrilateral elements.
1506 * The element has one gauss point where the result is calculated. In fact, the element has up to five gauss points
1507 * but GID does not support gauss points through the thickness of a Quadrilateral right now. Therefore, we will pick the middle gausspoint to show the results.
1508 * The point is placed on the natural position 0,0
1509 */
1510 out = new String(
1511 "GaussPoints \"Type" + type +
1512 "\" ElemType Quadrilateral \"MeshType" + type + "\"\n"
1513 );
1514 out += "Number Of Gauss Points: 1\n";
1515 out += "Nodes Not Included\n"; // There are no gauss points in the nodes.
1516 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.
1517 out += "0.0 0.0 \n";
1518 out += "End GaussPoints\n";
1519
1520 return out;
1521
1522 case RESULT_SUB_HEADER:
1523
1524 /* Print the subheader for the resultfile to initate each block of data from this element type
1525 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)
1526 Second param: Location of the data (1= on the nodes, 2= in the gauss points);
1527 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)
1528 Fourth param: Specify the name of the gauss point set that will be used "name"
1529 */
1530 out = new String(" 3 2 0 \"Type" + type + "\"\n");
1531
1532 return out;
1533
1534 case RESULT_STRESS_LOCAL:
1535
1536 /* Print the Gauss stresses for this element and the requested gauss point. */
1537 if (gpn == 0) {
1538 out = new String(number + " "); // Element number must start the first gauss point results
1539
1540 // The result is printed as [Sxx Syy Szz Sxy Syz Sxz]
1541 for (i = 0; i < 6; i++) {
1542 out += (stress[printed_integration_point].get(i, 0) + "\t");
1543 }
1544
1545 out += "\n";
1546 } else {
1547 out = new String("");
1548 }
1549
1550 return out; // Nothing more to print
1551
1552 case RESULT_STRAIN_LOCAL:
1553
1554 /* Print the Gauss strains for this element and the requested gauss point. */
1555 if (gpn == 0) {
1556 out = new String(number + " "); // Element number must start the first gauss point results
1557
1558 // The result is printed as [exx eyy ezz exy eyz exz]
1559 for (i = 0; i < 6; i++) {
1560 out += (strain[printed_integration_point].get(i, 0) + "\t");
1561 }
1562
1563 out += "\n";
1564 } else {
1565 out = new String("");
1566 }
1567
1568 return out; // Nothing more to print
1569
1570 default:
1571 return new String("");
1572 }
1573 }
1574
1575
1576 /**
1577 * This method is used to create the lines needed in the result file. The
1578 * method returns a string which is printed directly. However, due to the
1579 * fact that the line may be different depending on what is requested to
1580 * be printed and that the number of methods should be kept down, the
1581 * first parameter here is a control parameter. This parameter describes
1582 * what should be printed. The second parameter is a required input when
1583 * gauss point results are to be printed. Creation date: (09/12/01 %T)
1584 *
1585 * @param ctrl The control number to say if a header of result file is to
1586 * be printed.
1587 * @param gpn The gauss point number.
1588 *
1589 * @return java.lang.String
1590 */
print_Fembic(int ctrl, int gpn)1591 public String print_Fembic(int ctrl, int gpn) {
1592 String out;
1593
1594 switch (ctrl) {
1595
1596 case Element.MESH:
1597
1598 /* Print the element number and connected nodes */
1599 out = new String(
1600 number + "\t nodes = [" + node[0].getNumber() + ","
1601 + node[1].getNumber() + ","
1602 + node[2].getNumber() + ","
1603 + node[3].getNumber() + "]\t"
1604 + "T = " + thickness + "\t"
1605 + "material = " + material[0].getName() + "\t"
1606 );
1607
1608
1609 if (Factor_is_set)
1610 out += "factor = " + factor;
1611
1612 if (Friction_is_set)
1613 out += "friction = " + friction;
1614
1615 if (Contact == DISABLED)
1616 out += " contact = OFF";
1617
1618 if (Contact == EDGE)
1619 out += " contact = EDGE";
1620
1621 if (Contact == ADVANCED)
1622 out += " contact = ADVANCED";
1623
1624 if (Contact == ADVANCED_EDGE)
1625 out += " contact = ADVANCED_EDGE";
1626
1627 if (NIP_is_set)
1628 out += " nip = " + number_of_integration_points;
1629
1630 if (PIP_is_set)
1631 out += " pip = " + printed_integration_point;
1632
1633 if (shear_factor != 0)
1634 out += " shear_factor = " + shear_factor;
1635
1636 if (hourglass_control)
1637 out += " hourglass = ON";
1638
1639 if ( membrane_hourglass_coeff != 0)
1640 out += " mhc = " + membrane_hourglass_coeff;
1641
1642 if ( out_of_plane_hourglass_coeff != 0)
1643 out += " oophc = " + out_of_plane_hourglass_coeff;
1644
1645 if ( rotation_hourglass_coeff != 0)
1646 out += " rhc = " + rotation_hourglass_coeff;
1647
1648 if ( load != null)
1649 out += " load = " + load.getName();
1650
1651 if (!Thinning_is_enabled)
1652 out += " thinning = OFF";
1653
1654 out += "\n";
1655
1656 return out;
1657
1658
1659 default:
1660 return new String("");
1661 }
1662 }
1663
1664
1665 /**
1666 * Insert the method's description here. Creation date: (25/12/01 %T)
1667 */
setInitialConditions()1668 public void setInitialConditions()
1669 throws IllegalArgumentException
1670 {
1671 int i;
1672
1673 /* Make a local copy of the material object and keep it inside the element.
1674 * The reason for this is twofold.
1675 * 1. The material law is now tied to the element and will remember the element properties and history
1676 * 2. Since the material is cloned now, after Initialize, all the parameters are defined in the material object and are
1677 * now automatically transferred to the local copy.
1678 *
1679 * Note that one copy of the materiallaw is required for each intergration point in the element.
1680 */
1681 try {
1682 for (i = 0; i < number_of_integration_points; i++) {
1683 material[i] = (Material) material[0].copy();
1684 }
1685 } catch (CloneNotSupportedException e) {
1686 System.err.println("Object cannot clone");
1687 }
1688
1689 // Now call any necessary initialisations in the law.
1690 for (i = 0; i < number_of_integration_points; i++) {
1691 material[i].setInitialConditions();
1692 }
1693
1694 // Set up the index and weights of the integration points for a gaussian integration
1695 if (number_of_integration_points == 1) {
1696 z[0] = 0 * (thickness / 2.0);
1697 weight_factor[0] = 2.0;
1698 } else if (number_of_integration_points == 2) {
1699 z[0] = 0.577350269189626 * (thickness / 2.0);
1700 weight_factor[0] = 1.0;
1701 z[1] = -0.577350269189626 * (thickness / 2.0);
1702 weight_factor[1] = 1.0;
1703 } else if (number_of_integration_points == 3) {
1704 z[0] = 0.774596669241483 * (thickness / 2.0);
1705 weight_factor[0] = 5.0 / 9.0;
1706 z[1] = 0 * (thickness / 2.0);
1707 weight_factor[1] = 8.0 / 9.0;
1708 z[2] = -0.774596669241483 * (thickness / 2.0);
1709 weight_factor[2] = 5.0 / 9.0;
1710 } else if (number_of_integration_points == 4) {
1711 z[0] = 0.861136311594053 * (thickness / 2.0);
1712 weight_factor[0] = 0.347854845137454;
1713 z[1] = 0.339981043584856 * (thickness / 2.0);
1714 weight_factor[1] = 0.652145154862546;
1715 z[2] = -0.339981043584856 * (thickness / 2.0);
1716 weight_factor[2] = 0.652145154862546;
1717 z[3] = -0.861136311594053 * (thickness / 2.0);
1718 weight_factor[1] = 0.347854845137454;
1719 } else if (number_of_integration_points == 5) {
1720 z[0] = 0.90618 * (thickness / 2.0);
1721 weight_factor[0] = 0.236927;
1722 z[1] = 0.538469 * (thickness / 2.0);
1723 weight_factor[1] = 0.478629;
1724 z[2] = 0 * (thickness / 2.0);
1725 weight_factor[2] = 0.568889;
1726 z[3] = -0.538469 * (thickness / 2.0);
1727 weight_factor[3] = 0.478629;
1728 z[4] = -0.90618 * (thickness / 2.0);
1729 weight_factor[4] = 0.236927;
1730 }
1731
1732 // Define the hourglass vector
1733 hourglass_vector.set(0, 0, 1);
1734 hourglass_vector.set(1, 0, -1);
1735 hourglass_vector.set(2, 0, 1);
1736 hourglass_vector.set(3, 0, -1);
1737
1738 // Set initial conditions on the contact elements
1739 if ((Contact == BASIC) || (Contact == EDGE)) {
1740 internal_contact_element_1.setInitialConditions();
1741 internal_contact_element_2.setInitialConditions();
1742 }
1743
1744 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
1745 internal_contact_element_1.setInitialConditions();
1746 internal_contact_element_2.setInitialConditions();
1747 internal_contact_element_3.setInitialConditions();
1748 internal_contact_element_4.setInitialConditions();
1749
1750 }
1751
1752 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
1753 if (internal_contact_line_element_1 != null) {
1754 internal_contact_line_element_1.setInitialConditions();
1755 }
1756
1757 if (internal_contact_line_element_2 != null) {
1758 internal_contact_line_element_2.setInitialConditions();
1759 }
1760
1761 if (internal_contact_line_element_3 != null) {
1762 internal_contact_line_element_3.setInitialConditions();
1763 }
1764
1765 if (internal_contact_line_element_4 != null) {
1766 internal_contact_line_element_4.setInitialConditions();
1767 }
1768 }
1769
1770 // Print by default the middle integration point if nothing else specified
1771 if (! PIP_is_set) {
1772 printed_integration_point = number_of_integration_points / 2;
1773 }
1774 }
1775
1776 /**
1777 * This method calculates the local coordinate system for the element. The
1778 * method returns a handle to the system and this matrix is then stored in
1779 * the local_coordinate_system for later use. The matrix can be used in
1780 * the transformation of displacements and forces between the local and
1781 * global coordinates. However, in this element, the transformation is
1782 * made automatically in the matrix algebra of calculating the element
1783 * strains (they are derived in global directions) The element is assumed
1784 * to have the following node numbering y 3 I 2 I
1785 * I________ x 0 1
1786 */
updateLocalCoordinateSystem()1787 public void updateLocalCoordinateSystem() {
1788 this.calculateLocalBaseVectors(
1789 node[0].getX_pos(), node[0].getY_pos(), node[0].getZ_pos(),
1790 node[1].getX_pos(), node[1].getY_pos(), node[1].getZ_pos(),
1791 node[2].getX_pos(), node[2].getY_pos(), node[2].getZ_pos(),
1792 node[3].getX_pos(), node[3].getY_pos(), node[3].getZ_pos(),
1793 local_coordinate_system
1794 );
1795
1796 if ((Contact == BASIC) || (Contact == EDGE)) {
1797 internal_contact_element_1.updateLocalCoordinateSystem();
1798 internal_contact_element_2.updateLocalCoordinateSystem();
1799 }
1800
1801 if ((Contact == ADVANCED) || (Contact == ADVANCED_EDGE)) {
1802 internal_contact_element_1.updateLocalCoordinateSystem();
1803 internal_contact_element_2.updateLocalCoordinateSystem();
1804 internal_contact_element_3.updateLocalCoordinateSystem();
1805 internal_contact_element_4.updateLocalCoordinateSystem();
1806 }
1807
1808 if ((Contact == EDGE) || (Contact == ADVANCED_EDGE)) {
1809 if (internal_contact_line_element_1 != null) {
1810 internal_contact_line_element_1.updateLocalCoordinateSystem();
1811 }
1812
1813 if (internal_contact_line_element_2 != null) {
1814 internal_contact_line_element_2.updateLocalCoordinateSystem();
1815 }
1816
1817 if (internal_contact_line_element_3 != null) {
1818 internal_contact_line_element_3.updateLocalCoordinateSystem();
1819 }
1820
1821 if (internal_contact_line_element_4 != null) {
1822 internal_contact_line_element_4.updateLocalCoordinateSystem();
1823 }
1824 }
1825 }
1826
calculateLocalVariables()1827 private void calculateLocalVariables() {
1828 // Update the p matrix
1829 trash.set(0, 0, node[0].getX_pos());
1830 trash.set(1, 0, node[0].getY_pos());
1831 trash.set(2, 0, node[0].getZ_pos());
1832 trash.set(3, 0, node[1].getX_pos());
1833 trash.set(4, 0, node[1].getY_pos());
1834 trash.set(5, 0, node[1].getZ_pos());
1835 trash.set(6, 0, node[2].getX_pos());
1836 trash.set(7, 0, node[2].getY_pos());
1837 trash.set(8, 0, node[2].getZ_pos());
1838 trash.set(9, 0, node[3].getX_pos());
1839 trash.set(10, 0, node[3].getY_pos());
1840 trash.set(11, 0, node[3].getZ_pos());
1841
1842 // Transform into local coordinate system
1843 // p.setMatrix(0,2,0,0,local_coordinate_system.times(trash.getMatrix(0,2,0,0)));
1844 p.set(
1845 0, 0,
1846 (local_coordinate_system.get(0, 0) * trash.get(0, 0)) +
1847 (local_coordinate_system.get(0, 1) * trash.get(1, 0)) +
1848 (local_coordinate_system.get(0, 2) * trash.get(2, 0))
1849 );
1850 p.set(
1851 1, 0,
1852 (local_coordinate_system.get(1, 0) * trash.get(0, 0)) +
1853 (local_coordinate_system.get(1, 1) * trash.get(1, 0)) +
1854 (local_coordinate_system.get(1, 2) * trash.get(2, 0))
1855 );
1856 p.set(
1857 2, 0,
1858 (local_coordinate_system.get(2, 0) * trash.get(0, 0)) +
1859 (local_coordinate_system.get(2, 1) * trash.get(1, 0)) +
1860 (local_coordinate_system.get(2, 2) * trash.get(2, 0))
1861 );
1862
1863 // p.setMatrix(3,5,0,0,local_coordinate_system.times(trash.getMatrix(3,5,0,0)));
1864 p.set(
1865 3, 0,
1866 (local_coordinate_system.get(0, 0) * trash.get(3, 0)) +
1867 (local_coordinate_system.get(0, 1) * trash.get(4, 0)) +
1868 (local_coordinate_system.get(0, 2) * trash.get(5, 0))
1869 );
1870 p.set(
1871 4, 0,
1872 (local_coordinate_system.get(1, 0) * trash.get(3, 0)) +
1873 (local_coordinate_system.get(1, 1) * trash.get(4, 0)) +
1874 (local_coordinate_system.get(1, 2) * trash.get(5, 0))
1875 );
1876 p.set(
1877 5, 0,
1878 (local_coordinate_system.get(2, 0) * trash.get(3, 0)) +
1879 (local_coordinate_system.get(2, 1) * trash.get(4, 0)) +
1880 (local_coordinate_system.get(2, 2) * trash.get(5, 0))
1881 );
1882
1883 // p.setMatrix(6,8,0,0,local_coordinate_system.times(trash.getMatrix(6,8,0,0)));
1884 p.set(
1885 6, 0,
1886 (local_coordinate_system.get(0, 0) * trash.get(6, 0)) +
1887 (local_coordinate_system.get(0, 1) * trash.get(7, 0)) +
1888 (local_coordinate_system.get(0, 2) * trash.get(8, 0))
1889 );
1890 p.set(
1891 7, 0,
1892 (local_coordinate_system.get(1, 0) * trash.get(6, 0)) +
1893 (local_coordinate_system.get(1, 1) * trash.get(7, 0)) +
1894 (local_coordinate_system.get(1, 2) * trash.get(8, 0))
1895 );
1896 p.set(
1897 8, 0,
1898 (local_coordinate_system.get(2, 0) * trash.get(6, 0)) +
1899 (local_coordinate_system.get(2, 1) * trash.get(7, 0)) +
1900 (local_coordinate_system.get(2, 2) * trash.get(8, 0))
1901 );
1902
1903 // p.setMatrix(9,11,0,0,local_coordinate_system.times(trash.getMatrix(9,11,0,0)));
1904 p.set(
1905 9, 0,
1906 (local_coordinate_system.get(0, 0) * trash.get(9, 0)) +
1907 (local_coordinate_system.get(0, 1) * trash.get(10, 0)) +
1908 (local_coordinate_system.get(0, 2) * trash.get(11, 0))
1909 );
1910 p.set(
1911 10, 0,
1912 (local_coordinate_system.get(1, 0) * trash.get(9, 0)) +
1913 (local_coordinate_system.get(1, 1) * trash.get(10, 0)) +
1914 (local_coordinate_system.get(1, 2) * trash.get(11, 0))
1915 );
1916 p.set(
1917 11, 0,
1918 (local_coordinate_system.get(2, 0) * trash.get(9, 0)) +
1919 (local_coordinate_system.get(2, 1) * trash.get(10, 0)) +
1920 (local_coordinate_system.get(2, 2) * trash.get(11, 0))
1921 );
1922
1923 // Update the v matrix
1924 trash.set(0, 0, node[0].getX_vel());
1925 trash.set(1, 0, node[0].getY_vel());
1926 trash.set(2, 0, node[0].getZ_vel());
1927 trash.set(3, 0, node[1].getX_vel());
1928 trash.set(4, 0, node[1].getY_vel());
1929 trash.set(5, 0, node[1].getZ_vel());
1930 trash.set(6, 0, node[2].getX_vel());
1931 trash.set(7, 0, node[2].getY_vel());
1932 trash.set(8, 0, node[2].getZ_vel());
1933 trash.set(9, 0, node[3].getX_vel());
1934 trash.set(10, 0, node[3].getY_vel());
1935 trash.set(11, 0, node[3].getZ_vel());
1936
1937 // And transform into local coordinate system
1938 // v.setMatrix(0,2,0,0,local_coordinate_system.times(trash.getMatrix(0,2,0,0)));
1939 v.set(
1940 0, 0,
1941 (local_coordinate_system.get(0, 0) * trash.get(0, 0)) +
1942 (local_coordinate_system.get(0, 1) * trash.get(1, 0)) +
1943 (local_coordinate_system.get(0, 2) * trash.get(2, 0))
1944 );
1945 v.set(
1946 1, 0,
1947 (local_coordinate_system.get(1, 0) * trash.get(0, 0)) +
1948 (local_coordinate_system.get(1, 1) * trash.get(1, 0)) +
1949 (local_coordinate_system.get(1, 2) * trash.get(2, 0))
1950 );
1951 v.set(
1952 2, 0,
1953 (local_coordinate_system.get(2, 0) * trash.get(0, 0)) +
1954 (local_coordinate_system.get(2, 1) * trash.get(1, 0)) +
1955 (local_coordinate_system.get(2, 2) * trash.get(2, 0))
1956 );
1957
1958 // v.setMatrix(3,5,0,0,local_coordinate_system.times(trash.getMatrix(3,5,0,0)));
1959 v.set(
1960 3, 0,
1961 (local_coordinate_system.get(0, 0) * trash.get(3, 0)) +
1962 (local_coordinate_system.get(0, 1) * trash.get(4, 0)) +
1963 (local_coordinate_system.get(0, 2) * trash.get(5, 0))
1964 );
1965 v.set(
1966 4, 0,
1967 (local_coordinate_system.get(1, 0) * trash.get(3, 0)) +
1968 (local_coordinate_system.get(1, 1) * trash.get(4, 0)) +
1969 (local_coordinate_system.get(1, 2) * trash.get(5, 0))
1970 );
1971 v.set(
1972 5, 0,
1973 (local_coordinate_system.get(2, 0) * trash.get(3, 0)) +
1974 (local_coordinate_system.get(2, 1) * trash.get(4, 0)) +
1975 (local_coordinate_system.get(2, 2) * trash.get(5, 0))
1976 );
1977
1978 // v.setMatrix(6,8,0,0,local_coordinate_system.times(trash.getMatrix(6,8,0,0)));
1979 v.set(
1980 6, 0,
1981 (local_coordinate_system.get(0, 0) * trash.get(6, 0)) +
1982 (local_coordinate_system.get(0, 1) * trash.get(7, 0)) +
1983 (local_coordinate_system.get(0, 2) * trash.get(8, 0))
1984 );
1985 v.set(
1986 7, 0,
1987 (local_coordinate_system.get(1, 0) * trash.get(6, 0)) +
1988 (local_coordinate_system.get(1, 1) * trash.get(7, 0)) +
1989 (local_coordinate_system.get(1, 2) * trash.get(8, 0))
1990 );
1991 v.set(
1992 8, 0,
1993 (local_coordinate_system.get(2, 0) * trash.get(6, 0)) +
1994 (local_coordinate_system.get(2, 1) * trash.get(7, 0)) +
1995 (local_coordinate_system.get(2, 2) * trash.get(8, 0))
1996 );
1997
1998 // v.setMatrix(9,11,0,0,local_coordinate_system.times(trash.getMatrix(9,11,0,0)));
1999 v.set(
2000 9, 0,
2001 (local_coordinate_system.get(0, 0) * trash.get(9, 0)) +
2002 (local_coordinate_system.get(0, 1) * trash.get(10, 0)) +
2003 (local_coordinate_system.get(0, 2) * trash.get(11, 0))
2004 );
2005 v.set(
2006 10, 0,
2007 (local_coordinate_system.get(1, 0) * trash.get(9, 0)) +
2008 (local_coordinate_system.get(1, 1) * trash.get(10, 0)) +
2009 (local_coordinate_system.get(1, 2) * trash.get(11, 0))
2010 );
2011 v.set(
2012 11, 0,
2013 (local_coordinate_system.get(2, 0) * trash.get(9, 0)) +
2014 (local_coordinate_system.get(2, 1) * trash.get(10, 0)) +
2015 (local_coordinate_system.get(2, 2) * trash.get(11, 0))
2016 );
2017
2018 // Calculate the rotational velocity array
2019 trash.set(0, 0, node[0].getX_rot_vel());
2020 trash.set(1, 0, node[0].getY_rot_vel());
2021 trash.set(2, 0, node[0].getZ_rot_vel());
2022 trash.set(3, 0, node[1].getX_rot_vel());
2023 trash.set(4, 0, node[1].getY_rot_vel());
2024 trash.set(5, 0, node[1].getZ_rot_vel());
2025 trash.set(6, 0, node[2].getX_rot_vel());
2026 trash.set(7, 0, node[2].getY_rot_vel());
2027 trash.set(8, 0, node[2].getZ_rot_vel());
2028 trash.set(9, 0, node[3].getX_rot_vel());
2029 trash.set(10, 0, node[3].getY_rot_vel());
2030 trash.set(11, 0, node[3].getZ_rot_vel());
2031
2032 // And transform into local coordinate system
2033 // a.setMatrix(0,2,0,0,local_coordinate_system.times(trash.getMatrix(0,2,0,0)));
2034 a.set(
2035 0, 0,
2036 (local_coordinate_system.get(0, 0) * trash.get(0, 0)) +
2037 (local_coordinate_system.get(0, 1) * trash.get(1, 0)) +
2038 (local_coordinate_system.get(0, 2) * trash.get(2, 0))
2039 );
2040 a.set(
2041 1, 0,
2042 (local_coordinate_system.get(1, 0) * trash.get(0, 0)) +
2043 (local_coordinate_system.get(1, 1) * trash.get(1, 0)) +
2044 (local_coordinate_system.get(1, 2) * trash.get(2, 0))
2045 );
2046 a.set(
2047 2, 0,
2048 (local_coordinate_system.get(2, 0) * trash.get(0, 0)) +
2049 (local_coordinate_system.get(2, 1) * trash.get(1, 0)) +
2050 (local_coordinate_system.get(2, 2) * trash.get(2, 0))
2051 );
2052
2053 // a.setMatrix(3,5,0,0,local_coordinate_system.times(trash.getMatrix(3,5,0,0)));
2054 a.set(
2055 3, 0,
2056 (local_coordinate_system.get(0, 0) * trash.get(3, 0)) +
2057 (local_coordinate_system.get(0, 1) * trash.get(4, 0)) +
2058 (local_coordinate_system.get(0, 2) * trash.get(5, 0))
2059 );
2060 a.set(
2061 4, 0,
2062 (local_coordinate_system.get(1, 0) * trash.get(3, 0)) +
2063 (local_coordinate_system.get(1, 1) * trash.get(4, 0)) +
2064 (local_coordinate_system.get(1, 2) * trash.get(5, 0))
2065 );
2066 a.set(
2067 5, 0,
2068 (local_coordinate_system.get(2, 0) * trash.get(3, 0)) +
2069 (local_coordinate_system.get(2, 1) * trash.get(4, 0)) +
2070 (local_coordinate_system.get(2, 2) * trash.get(5, 0))
2071 );
2072
2073 // a.setMatrix(6,8,0,0,local_coordinate_system.times(trash.getMatrix(6,8,0,0)));
2074 a.set(
2075 6, 0,
2076 (local_coordinate_system.get(0, 0) * trash.get(6, 0)) +
2077 (local_coordinate_system.get(0, 1) * trash.get(7, 0)) +
2078 (local_coordinate_system.get(0, 2) * trash.get(8, 0))
2079 );
2080 a.set(
2081 7, 0,
2082 (local_coordinate_system.get(1, 0) * trash.get(6, 0)) +
2083 (local_coordinate_system.get(1, 1) * trash.get(7, 0)) +
2084 (local_coordinate_system.get(1, 2) * trash.get(8, 0))
2085 );
2086 a.set(
2087 8, 0,
2088 (local_coordinate_system.get(2, 0) * trash.get(6, 0)) +
2089 (local_coordinate_system.get(2, 1) * trash.get(7, 0)) +
2090 (local_coordinate_system.get(2, 2) * trash.get(8, 0))
2091 );
2092
2093 // a.setMatrix(9,11,0,0,local_coordinate_system.times(trash.getMatrix(9,11,0,0)));
2094 a.set(
2095 9, 0,
2096 (local_coordinate_system.get(0, 0) * trash.get(9, 0)) +
2097 (local_coordinate_system.get(0, 1) * trash.get(10, 0)) +
2098 (local_coordinate_system.get(0, 2) * trash.get(11, 0))
2099 );
2100 a.set(
2101 10, 0,
2102 (local_coordinate_system.get(1, 0) * trash.get(9, 0)) +
2103 (local_coordinate_system.get(1, 1) * trash.get(10, 0)) +
2104 (local_coordinate_system.get(1, 2) * trash.get(11, 0))
2105 );
2106 a.set(
2107 11, 0,
2108 (local_coordinate_system.get(2, 0) * trash.get(9, 0)) +
2109 (local_coordinate_system.get(2, 1) * trash.get(10, 0)) +
2110 (local_coordinate_system.get(2, 2) * trash.get(11, 0))
2111 );
2112
2113 /* Calculate element area
2114 * Normal formula works for small distorsions: area = (1/2)*(X31*Y42 + X24*Y31)
2115 * Note that for large deformations, the area will become very small as the element
2116 * reaches self contact. The area is then corrected and recalculated in the
2117 * calculateStrain() method if the ADVANCED contact option is used.
2118 */
2119 area = 0.5 * (((p.get(6, 0) - p.get(0, 0)) * (p.get(10, 0) -
2120 p.get(4, 0))) +
2121 ((p.get(3, 0) - p.get(9, 0)) * (p.get(7, 0) - p.get(1, 0))));
2122
2123 }
2124
2125
2126
2127 /**
2128 * This method is called by any internal nodes that this element uses.
2129 * This is the case if advanced or advanced_edge contact options have been
2130 * used. The method is called by the internal node when the calculateNewPosition
2131 * is run for all the nodes in the main loop.
2132 */
2133
setInternalNodePosition()2134 public void setInternalNodePosition() {
2135
2136 // Update the position of the local node if needed
2137 middle_node.setX_pos_orig(
2138 (node[0].getX_pos() + node[1].getX_pos() + node[2].getX_pos() +
2139 node[3].getX_pos()) / 4
2140 );
2141 middle_node.setY_pos_orig(
2142 (node[0].getY_pos() + node[1].getY_pos() + node[2].getY_pos() +
2143 node[3].getY_pos()) / 4
2144 );
2145 middle_node.setZ_pos_orig(
2146 (node[0].getZ_pos() + node[1].getZ_pos() + node[2].getZ_pos() +
2147 node[3].getZ_pos()) / 4
2148 );
2149 }
2150
2151
getInternalNode()2152 public Node getInternalNode() {
2153 return middle_node;
2154 }
2155
2156 /**
2157 * Is used to deactivate the element when fractured.
2158 */
deActivate()2159 public void deActivate() {
2160 super.deActivate();
2161 if (internal_contact_element_1 != null) internal_contact_element_1.deActivate();
2162 if (internal_contact_element_2 != null) internal_contact_element_2.deActivate();
2163 if (internal_contact_element_3 != null) internal_contact_element_3.deActivate();
2164 if (internal_contact_element_4 != null) internal_contact_element_4.deActivate();
2165 if (internal_contact_line_element_1 != null) internal_contact_line_element_1.deActivate();
2166 if (internal_contact_line_element_2 != null) internal_contact_line_element_2.deActivate();
2167 if (internal_contact_line_element_3 != null) internal_contact_line_element_3.deActivate();
2168 if (internal_contact_line_element_4 != null) internal_contact_line_element_4.deActivate();
2169 if (middle_node != null) middle_node.deActivate();
2170 }
2171
2172
2173 }
2174
2175