00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00024
00025
00026 #ifndef TWENTYSEVENNODEBRICK_CPP
00027 #define TWENTYSEVENNODEBRICK_CPP
00028
00029 #include <NDMaterial.h>
00030 #include <Matrix.h>
00031 #include <Vector.h>
00032 #include <ID.h>
00033 #include <Renderer.h>
00034 #include <Domain.h>
00035 #include <string.h>
00036 #include <Information.h>
00037 #include <Channel.h>
00038 #include <FEM_ObjectBroker.h>
00039 #include <ElementResponse.h>
00040
00041 #include <TwentySevenNodeBrick.h>
00042 #include <ElementalLoad.h>
00043 #define FixedOrder 3
00044
00045
00046 Matrix TwentySevenNodeBrick::K(81, 81);
00047 Matrix TwentySevenNodeBrick::C(81, 81);
00048 Matrix TwentySevenNodeBrick::M(81, 81);
00049 Vector TwentySevenNodeBrick::P(81);
00050 Vector InfoMoment(109+3);
00051 Vector InfoPlastic(FixedOrder*FixedOrder*FixedOrder*4+1);
00052 Vector InfoStress(FixedOrder*FixedOrder*FixedOrder*6+1);
00053 Vector GaussCoord(FixedOrder*FixedOrder*FixedOrder*3+1);
00054 Vector Info_pq2(2);
00055
00056
00057
00058
00059 TwentySevenNodeBrick::TwentySevenNodeBrick(int element_number,
00060 int node_numb_1, int node_numb_2, int node_numb_3, int node_numb_4,
00061 int node_numb_5, int node_numb_6, int node_numb_7, int node_numb_8,
00062 int node_numb_9, int node_numb_10, int node_numb_11, int node_numb_12,
00063 int node_numb_13, int node_numb_14, int node_numb_15, int node_numb_16,
00064 int node_numb_17, int node_numb_18, int node_numb_19, int node_numb_20,
00065 int node_numb_21, int node_numb_22, int node_numb_23, int node_numb_24,
00066 int node_numb_25, int node_numb_26, int node_numb_27,
00067 NDMaterial * Globalmmodel, double b1, double b2,double b3,
00068 double r, double p)
00069
00070 :Element(element_number, ELE_TAG_TwentySevenNodeBrick ),
00071 connectedExternalNodes(27), Ki(0), Q(81), bf(3),
00072 rho(r), pressure(p)
00073 {
00074
00075 bf(0) = b1;
00076 bf(1) = b2;
00077 bf(2) = b3;
00078
00079 determinant_of_Jacobian = 0.0;
00080
00081
00082
00083
00084 r_integration_order = FixedOrder;
00085 s_integration_order = FixedOrder;
00086 t_integration_order = FixedOrder;
00087
00088
00089
00090
00091 int total_number_of_Gauss_points = r_integration_order*s_integration_order*t_integration_order;
00092
00093
00094 if ( total_number_of_Gauss_points != 0 )
00095 {
00096 matpoint = new MatPoint3D * [total_number_of_Gauss_points];
00097 }
00098 else
00099 {
00100 matpoint = 0;
00101 }
00103 short where = 0;
00104
00105 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
00106 {
00107 double r = get_Gauss_p_c( r_integration_order, GP_c_r );
00108 double rw = get_Gauss_p_w( r_integration_order, GP_c_r );
00109
00110 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
00111 {
00112 double s = get_Gauss_p_c( s_integration_order, GP_c_s );
00113 double sw = get_Gauss_p_w( s_integration_order, GP_c_s );
00114
00115 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
00116 {
00117 double t = get_Gauss_p_c( t_integration_order, GP_c_t );
00118 double tw = get_Gauss_p_w( t_integration_order, GP_c_t );
00119
00120
00121
00122 where =
00123 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 matpoint[where] = new MatPoint3D(GP_c_r,
00141 GP_c_s,
00142 GP_c_t,
00143 r, s, t,
00144 rw, sw, tw,
00145
00146 Globalmmodel);
00147
00148
00149
00150
00151 }
00152 }
00153 }
00154
00155
00156 connectedExternalNodes( 0) = node_numb_1;
00157 connectedExternalNodes( 1) = node_numb_2;
00158 connectedExternalNodes( 2) = node_numb_3;
00159 connectedExternalNodes( 3) = node_numb_4;
00160 connectedExternalNodes( 4) = node_numb_5;
00161 connectedExternalNodes( 5) = node_numb_6;
00162 connectedExternalNodes( 6) = node_numb_7;
00163 connectedExternalNodes( 7) = node_numb_8;
00164
00165 connectedExternalNodes( 8) = node_numb_9;
00166 connectedExternalNodes( 9) = node_numb_10;
00167 connectedExternalNodes(10) = node_numb_11;
00168 connectedExternalNodes(11) = node_numb_12;
00169
00170 connectedExternalNodes(12) = node_numb_13;
00171 connectedExternalNodes(13) = node_numb_14;
00172 connectedExternalNodes(14) = node_numb_15;
00173 connectedExternalNodes(15) = node_numb_16;
00174
00175 connectedExternalNodes(16) = node_numb_17;
00176 connectedExternalNodes(17) = node_numb_18;
00177 connectedExternalNodes(18) = node_numb_19;
00178 connectedExternalNodes(19) = node_numb_20;
00179
00180 connectedExternalNodes(20) = node_numb_21;
00181 connectedExternalNodes(21) = node_numb_22;
00182 connectedExternalNodes(22) = node_numb_23;
00183 connectedExternalNodes(23) = node_numb_24;
00184 connectedExternalNodes(24) = node_numb_25;
00185 connectedExternalNodes(25) = node_numb_26;
00186 connectedExternalNodes(26) = node_numb_27;
00187
00188 for (int i=0; i<27; i++)
00189 theNodes[i] = 0;
00190
00191 nodes_in_brick = 27;
00192
00193 }
00194
00195
00196 TwentySevenNodeBrick::TwentySevenNodeBrick ():Element(0, ELE_TAG_TwentySevenNodeBrick ),
00197 connectedExternalNodes(27), Ki(0), Q(81), bf(3), rho(0.0), pressure(0.0), mmodel(0)
00198 {
00199 matpoint = 0;
00200
00201 for (int i=0; i<27; i++)
00202 theNodes[i] = 0;
00203 }
00204
00205
00206
00207
00208
00211 TwentySevenNodeBrick::~TwentySevenNodeBrick ()
00212 {
00213
00214 int total_number_of_Gauss_points = r_integration_order*s_integration_order*t_integration_order;
00215
00216 for (int i = 0; i < total_number_of_Gauss_points; i++)
00217 {
00218
00219 if (matpoint[i])
00220 delete matpoint[i];
00221 }
00222
00223
00224 if (matpoint)
00225 delete [] matpoint;
00226
00227 if (Ki != 0)
00228 delete Ki;
00229
00230 }
00231
00235 void TwentySevenNodeBrick::incremental_Update()
00236 {
00237 double r = 0.0;
00238
00239 double s = 0.0;
00240
00241 double t = 0.0;
00242
00243
00244 short where = 0;
00245
00246
00247
00248
00249 int dh_dim[] = {27,3};
00250 tensor dh(2, dh_dim, 0.0);
00251
00252
00253 static int disp_dim[] = {27,3};
00254 tensor incremental_displacements(2,disp_dim,0.0);
00255
00256 straintensor incremental_strain;
00257
00258 tensor Jacobian;
00259 tensor JacobianINV;
00260 tensor dhGlobal;
00261
00262 incremental_displacements = incr_disp();
00263
00264 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
00265 {
00266 r = get_Gauss_p_c( r_integration_order, GP_c_r );
00267
00268 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
00269 {
00270 s = get_Gauss_p_c( s_integration_order, GP_c_s );
00271
00272 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
00273 {
00274 t = get_Gauss_p_c( t_integration_order, GP_c_t );
00275
00276
00277
00278 where =
00279 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
00280
00281 dh = dh_drst_at(r,s,t);
00282
00283 Jacobian = Jacobian_3D(dh);
00284
00285
00286 JacobianINV = Jacobian_3Dinv(dh);
00287
00288
00289
00290
00291
00292
00293 dhGlobal = dh("ij") * JacobianINV("kj");
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 incremental_strain =
00308 (dhGlobal("ib")*incremental_displacements("ia")).symmetrize11();
00309 incremental_strain.null_indices();
00310
00311
00312
00313
00314
00315
00316 if ( ! ( (matpoint[where]->matmodel)->setTrialStrainIncr( incremental_strain)) )
00317 opserr << "TwentySevenNodeBrick::incremental_Update (tag: " << this->getTag() << "), not converged\n";
00318
00319 }
00320 }
00321 }
00322 }
00323
00324
00325
00326
00327
00328 tensor TwentySevenNodeBrick::H_3D(double r1, double r2, double r3)
00329 {
00330
00331 int dimension[] = {81,3};
00332
00333 tensor H(2, dimension, 0.0);
00334
00335
00336 H.val(1,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00337 H.val(2,2)=H.val(1,1);
00338 H.val(3,3)=H.val(1,1);
00339
00340
00341 H.val(4,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00342 H.val(5,2)=H.val(4,1);
00343 H.val(6,3)=H.val(4,1);
00344
00345
00346 H.val(7,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00347 H.val(8,2)=H.val(7,1);
00348 H.val(9,3)=H.val(7,1);
00349
00350
00351 H.val(10,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00352 H.val(11,2)=H.val(10,1);
00353 H.val(12,3)=H.val(10,1);
00354
00355
00356 H.val(13,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00357 H.val(14,2)=H.val(13,1);
00358 H.val(15,3)=H.val(13,1);
00359
00360
00361 H.val(16,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00362 H.val(17,2)=H.val(16,1);
00363 H.val(18,3)=H.val(16,1);
00364
00365
00366 H.val(19,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00367 H.val(20,2)=H.val(19,1);
00368 H.val(21,3)=H.val(19,1);
00369
00370
00371
00372 H.val(22,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00373 H.val(23,2)=H.val(22,1);
00374 H.val(24,3)=H.val(22,1);
00375
00376
00377
00378 H.val(25,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00379 H.val(26,2)=H.val(25,1);
00380 H.val(27,3)=H.val(25,1);
00381
00382
00383
00384 H.val(28,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00385 H.val(29,2)=H.val(28,1);
00386 H.val(30,3)=H.val(28,1);
00387
00388
00389
00390 H.val(31,1)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00391 H.val(32,2)=H.val(31,1);
00392 H.val(33,3)=H.val(31,1);
00393
00394
00395
00396 H.val(34,1)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00397 H.val(35,2)=H.val(34,1);
00398 H.val(36,3)=H.val(34,1);
00399
00400
00401
00402 H.val(37,1)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00403 H.val(38,2)=H.val(37,1);
00404 H.val(39,3)=H.val(37,1);
00405
00406
00407
00408 H.val(40,1)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00409 H.val(41,2)=H.val(40,1);
00410 H.val(42,3)=H.val(40,1);
00411
00412
00413
00414 H.val(43,1)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00415 H.val(44,2)=H.val(43,1);
00416 H.val(45,3)=H.val(43,1);
00417
00418
00419
00420 H.val(46,1)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00421 H.val(47,2)=H.val(46,1);
00422 H.val(48,3)=H.val(46,1);
00423
00424
00425
00426 H.val(49,1)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00427 H.val(50,2)=H.val(49,1);
00428 H.val(51,3)=H.val(49,1);
00429
00430
00431
00432 H.val(52,1)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00433 H.val(53,2)=H.val(52,1);
00434 H.val(54,3)=H.val(52,1);
00435
00436
00437
00438 H.val(55,1)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00439 H.val(56,2)=H.val(55,1);
00440 H.val(57,3)=H.val(55,1);
00441
00442
00443
00444 H.val(58,1)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00445 H.val(59,2)=H.val(58,1);
00446 H.val(60,3)=H.val(58,1);
00447
00448
00449
00450 H.val(61,1)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00451 H.val(62,2)=H.val(61,1);
00452 H.val(63,3)=H.val(61,1);
00453
00454
00455
00456 H.val(64,1)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*(1.0-r3*r3);
00457 H.val(65,2)=H.val(64,1);
00458 H.val(66,3)=H.val(64,1);
00459
00460
00461
00462 H.val(67,1)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00463 H.val(68,2)=H.val(67,1);
00464 H.val(69,3)=H.val(67,1);
00465
00466
00467
00468 H.val(70,1)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*(1.0-r3*r3);
00469 H.val(71,2)=H.val(70,1);
00470 H.val(72,3)=H.val(70,1);
00471
00472
00473
00474 H.val(73,1)=(1.0-r1*r1)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00475 H.val(74,2)=H.val(73,1);
00476 H.val(75,3)=H.val(73,1);
00477
00478
00479
00480 H.val(76,1)=(1.0-r1*r1)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00481 H.val(77,2)=H.val(76,1);
00482 H.val(78,3)=H.val(76,1);
00483
00484
00485
00486 H.val(79,1)=(1.0-r1*r1)*(1.0-r2*r2)*(1.0-r3*r3);
00487 H.val(80,2)=H.val(79,1);
00488 H.val(81,3)=H.val(79,1);
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 return H;
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664 tensor TwentySevenNodeBrick::dh_drst_at(double r1, double r2, double r3)
00665 {
00666
00667 int dimensions[] = {27,3};
00668 tensor dh(2, dimensions, 0.0);
00669
00670
00671
00672 dh.val(1,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00673 dh.val(1,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2+1.0)*0.5*r3*(r3+1.0);
00674 dh.val(1,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*0.5*(2.0*r3+1.0);
00675
00676
00677 dh.val(2,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00678 dh.val(2,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2+1.0)*0.5*r3*(r3+1.0);
00679 dh.val(2,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*0.5*(2.0*r3+1.0);
00680
00681
00682 dh.val(3,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00683 dh.val(3,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2-1.0)*0.5*r3*(r3+1.0);
00684 dh.val(3,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*0.5*(2.0*r3+1.0);
00685
00686
00687 dh.val(4,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00688 dh.val(4,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2-1.0)*0.5*r3*(r3+1.0);
00689 dh.val(4,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*0.5*(2.0*r3+1.0);
00690
00691
00692 dh.val(5,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00693 dh.val(5,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2+1.0)*0.5*r3*(r3-1.0);
00694 dh.val(5,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*0.5*(2.0*r3-1.0);
00695
00696
00697 dh.val(6,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00698 dh.val(6,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2+1.0)*0.5*r3*(r3-1.0);
00699 dh.val(6,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*0.5*(2.0*r3-1.0);
00700
00701
00702 dh.val(7,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00703 dh.val(7,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2-1.0)*0.5*r3*(r3-1.0);
00704 dh.val(7,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*0.5*(2.0*r3-1.0);
00705
00706
00707
00708 dh.val(8,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00709 dh.val(8,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2-1.0)*0.5*r3*(r3-1.0);
00710 dh.val(8,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*0.5*(2.0*r3-1.0);
00711
00712
00713
00714 dh.val(9,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00715 dh.val(9,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2+1.0)*(1.0-r3*r3);
00716 dh.val(9,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2+1.0)*(-2.0*r3);
00717
00718
00719
00720 dh.val(10,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00721 dh.val(10,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2+1.0)*(1.0-r3*r3);
00722 dh.val(10,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2+1.0)*(-2.0*r3);
00723
00724
00725
00726 dh.val(11,1)=0.5*(2.0*r1-1.0)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00727 dh.val(11,2)=0.5*r1*(r1-1.0)*0.5*(2.0*r2-1.0)*(1.0-r3*r3);
00728 dh.val(11,3)=0.5*r1*(r1-1.0)*0.5*r2*(r2-1.0)*(-2.0*r3);
00729
00730
00731
00732 dh.val(12,1)=0.5*(2.0*r1+1.0)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00733 dh.val(12,2)=0.5*r1*(r1+1.0)*0.5*(2.0*r2-1.0)*(1.0-r3*r3);
00734 dh.val(12,3)=0.5*r1*(r1+1.0)*0.5*r2*(r2-1.0)*(-2.0*r3);
00735
00736
00737
00738 dh.val(13,1)=(-2.0*r1)*0.5*r2*(r2+1.0)*0.5*r3*(r3+1.0);
00739 dh.val(13,2)=(1.0-r1*r1)*0.5*(2.0*r2+1.0)*0.5*r3*(r3+1.0);
00740 dh.val(13,3)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*0.5*(2.0*r3+1.0);
00741
00742
00743
00744 dh.val(14,1)=0.5*(2.0*r1-1.0)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00745 dh.val(14,2)=0.5*r1*(r1-1.0)*(-2.0*r2)*0.5*r3*(r3+1.0);
00746 dh.val(14,3)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*0.5*(2.0*r3+1.0);
00747
00748
00749
00750 dh.val(15,1)=(-2.0*r1)*0.5*r2*(r2-1.0)*0.5*r3*(r3+1.0);
00751 dh.val(15,2)=(1.0-r1*r1)*0.5*(2.0*r2-1.0)*0.5*r3*(r3+1.0);
00752 dh.val(15,3)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*0.5*(2.0*r3+1.0);
00753
00754
00755
00756 dh.val(16,1)=0.5*(2.0*r1+1.0)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00757 dh.val(16,2)=0.5*r1*(r1+1.0)*(-2.0*r2)*0.5*r3*(r3+1.0);
00758 dh.val(16,3)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*0.5*(2.0*r3+1.0);
00759
00760
00761
00762 dh.val(17,1)=(-2.0*r1)*0.5*r2*(r2+1.0)*0.5*r3*(r3-1.0);
00763 dh.val(17,2)=(1.0-r1*r1)*0.5*(2.0*r2+1.0)*0.5*r3*(r3-1.0);
00764 dh.val(17,3)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*0.5*(2.0*r3-1.0);
00765
00766
00767
00768 dh.val(18,1)=0.5*(2.0*r1-1.0)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00769 dh.val(18,2)=0.5*r1*(r1-1.0)*(-2.0*r2)*0.5*r3*(r3-1.0);
00770 dh.val(18,3)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*0.5*(2.0*r3-1.0);
00771
00772
00773
00774 dh.val(19,1)=(-2.0*r1)*0.5*r2*(r2-1.0)*0.5*r3*(r3-1.0);
00775 dh.val(19,2)=(1.0-r1*r1)*0.5*(2.0*r2-1.0)*0.5*r3*(r3-1.0);
00776 dh.val(19,3)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*0.5*(2.0*r3-1.0);
00777
00778
00779
00780 dh.val(20,1)=0.5*(2.0*r1+1.0)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00781 dh.val(20,2)=0.5*r1*(r1+1.0)*(-2.0*r2)*0.5*r3*(r3-1.0);
00782 dh.val(20,3)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*0.5*(2.0*r3-1.0);
00783
00784
00785
00786 dh.val(21,1)=(-2.0*r1)*0.5*r2*(r2+1.0)*(1.0-r3*r3);
00787 dh.val(21,2)=(1.0-r1*r1)*0.5*(2.0*r2+1.0)*(1.0-r3*r3);
00788 dh.val(21,3)=(1.0-r1*r1)*0.5*r2*(r2+1.0)*(-2.0*r3);
00789
00790
00791
00792 dh.val(22,1)=0.5*(2.0*r1-1.0)*(1.0-r2*r2)*(1.0-r3*r3);
00793 dh.val(22,2)=0.5*r1*(r1-1.0)*(-2.0*r2)*(1.0-r3*r3);
00794 dh.val(22,3)=0.5*r1*(r1-1.0)*(1.0-r2*r2)*(-2.0*r3);
00795
00796
00797
00798 dh.val(23,1)=(-2.0*r1)*0.5*r2*(r2-1.0)*(1.0-r3*r3);
00799 dh.val(23,2)=(1.0-r1*r1)*0.5*(2.0*r2-1.0)*(1.0-r3*r3);
00800 dh.val(23,3)=(1.0-r1*r1)*0.5*r2*(r2-1.0)*(-2.0*r3);
00801
00802
00803
00804 dh.val(24,1)=0.5*(2.0*r1+1.0)*(1.0-r2*r2)*(1.0-r3*r3);
00805 dh.val(24,2)=0.5*r1*(r1+1.0)*(-2.0*r2)*(1.0-r3*r3);
00806 dh.val(24,3)=0.5*r1*(r1+1.0)*(1.0-r2*r2)*(-2.0*r3);
00807
00808
00809
00810 dh.val(25,1)=(-2.0*r1)*(1.0-r2*r2)*0.5*r3*(r3+1.0);
00811 dh.val(25,2)=(1.0-r1*r1)*(-2.0*r2)*0.5*r3*(r3+1.0);
00812 dh.val(25,3)=(1.0-r1*r1)*(1.0-r2*r2)*0.5*(2.0*r3+1.0);
00813
00814
00815
00816 dh.val(26,1)=(-2.0*r1)*(1.0-r2*r2)*0.5*r3*(r3-1.0);
00817 dh.val(26,2)=(1.0-r1*r1)*(-2.0*r2)*0.5*r3*(r3-1.0);
00818 dh.val(26,3)=(1.0-r1*r1)*(1.0-r2*r2)*0.5*(2.0*r3-1.0);
00819
00820
00821
00822 dh.val(27,1)=(-2.0*r1)*(1.0-r2*r2)*(1.0-r3*r3);
00823 dh.val(27,2)=(1.0-r1*r1)*(-2.0*r2)*(1.0-r3*r3);
00824 dh.val(27,3)=(1.0-r1*r1)*(1.0-r2*r2)*(-2.0*r3);
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914 return dh;
00915 }
00916
00917
00918
00920 TwentySevenNodeBrick & TwentySevenNodeBrick::operator[](int subscript)
00921 {
00922 return ( *(this+subscript) );
00923 }
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00940 tensor TwentySevenNodeBrick::getStiffnessTensor(void)
00941 {
00942 int K_dim[] = {27,3,3,27};
00943 tensor Kk(4,K_dim,0.0);
00944 tensor Kkt(4,K_dim,0.0);
00945
00946
00947
00948
00949 double r = 0.0;
00950 double rw = 0.0;
00951 double s = 0.0;
00952 double sw = 0.0;
00953 double t = 0.0;
00954 double tw = 0.0;
00955
00956 short where = 0;
00957 double weight = 0.0;
00958
00959 int dh_dim[] = {27,3};
00960 tensor dh(2, dh_dim, 0.0);
00961
00962
00963 tensor Constitutive;
00964
00965 double det_of_Jacobian = 0.0;
00966
00967 static int disp_dim[] = {27,3};
00968 tensor incremental_displacements(2,disp_dim,0.0);
00969
00970 straintensor incremental_strain;
00971
00972
00973 tensor Jacobian;
00974 tensor JacobianINV;
00975 tensor JacobianINVtemp;
00976 tensor dhGlobal;
00977
00978 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
00979 {
00980 r = get_Gauss_p_c( r_integration_order, GP_c_r );
00981 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
00982 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
00983 {
00984 s = get_Gauss_p_c( s_integration_order, GP_c_s );
00985 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
00986 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
00987 {
00988 t = get_Gauss_p_c( t_integration_order, GP_c_t );
00989 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
00990
00991
00992 where =
00993 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
00994
00995 dh = dh_drst_at(r,s,t);
00996
00997
00998 Jacobian = Jacobian_3D(dh);
00999
01000 JacobianINV = Jacobian_3Dinv(dh);
01001
01002
01003 det_of_Jacobian = Jacobian.determinant();
01004
01006 dhGlobal = dh("ij") * JacobianINV("kj");
01007
01008
01009
01010 weight = rw * sw * tw * det_of_Jacobian;
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089 Constitutive = (matpoint[where]->matmodel)->getTangentTensor();
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116 Kkt = dhGlobal("ib")*Constitutive("abcd");
01117 Kk = Kk + Kkt("aicd")*dhGlobal("jd")*weight;
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129 }
01130 }
01131 }
01132
01133
01134 return Kk;
01135 }
01136
01137
01138
01139
01140 void TwentySevenNodeBrick::set_strain_stress_tensor(FILE *fp, double * u)
01141 {
01142 int dh_dim[] = {27,3};
01143 tensor dh(2, dh_dim, 0.0);
01144
01145
01146 tensor Constitutive;
01147 double r = 0.0;
01148 double s = 0.0;
01149 double t = 0.0;
01150 int where = 0;
01151
01152 double det_of_Jacobian;
01153
01154 straintensor strain;
01155 stresstensor stress;
01156
01157 tensor Jacobian;
01158 tensor JacobianINV;
01159 tensor dhGlobal;
01160
01161
01162 static int disp_dim[] = {27,3};
01163 tensor total_displacements(2,disp_dim,0.0);
01164
01165 total_displacements = total_disp(fp, u);
01166
01167 ::printf("\n displacement(x-y-z) at GAUSS pt %d \n\n", where+1);
01168 for (int ii=1; ii<=27;ii++)
01169 {
01170 ::printf("Global# %d Local#%d %+0.5e %+0.5e %+0.5e\n",
01171
01172 connectedExternalNodes(ii-1),
01173 ii,total_displacements.val(ii,1),
01174 total_displacements.val(ii,2),
01175 total_displacements.val(ii,3));
01176 }
01177
01178 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
01179 {
01180 r = get_Gauss_p_c( r_integration_order, GP_c_r );
01181 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
01182 {
01183 s = get_Gauss_p_c( s_integration_order, GP_c_s );
01184 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
01185 {
01186 t = get_Gauss_p_c( t_integration_order, GP_c_t );
01187
01188
01189 where =
01190 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
01191
01192 dh = dh_drst_at(r,s,t);
01193
01194 Jacobian = Jacobian_3D(dh);
01195
01196 JacobianINV = Jacobian_3Dinv(dh);
01197
01198 det_of_Jacobian = Jacobian.determinant();
01199
01200
01201 dhGlobal = dh("ij") * JacobianINV("kj");
01202
01203
01204 strain = (dhGlobal("ib")*total_displacements("ia")).symmetrize11();
01205 strain.null_indices();
01206
01207
01208
01209
01210
01211
01212 Constitutive = ( matpoint[where]->matmodel)->getTangentTensor();
01213
01214 stress = Constitutive("ijkl") * strain("kl");
01215 stress.null_indices();
01216
01217 ::printf("\n strain tensor at GAUSS point %d \n", where+1);
01218 strain.reportshort("");
01219 ::printf("\n stress tensor at GAUSS point %d \n", where+1);
01220 stress.reportshort("");
01221
01222
01223 }
01224 }
01225 }
01226 }
01227
01228
01230
01231 tensor TwentySevenNodeBrick::getMassTensor(void)
01232 {
01233
01234 int M_dim[] = {81,81};
01235 tensor Mm(2,M_dim,0.0);
01236
01237 double r = 0.0;
01238 double rw = 0.0;
01239 double s = 0.0;
01240 double sw = 0.0;
01241 double t = 0.0;
01242 double tw = 0.0;
01243
01244 short where = 0;
01245 double weight = 0.0;
01246
01247 int dh_dim[] = {27,3};
01248
01249 tensor dh(2, dh_dim, 0.0);
01250
01251 int h_dim[] = {81,3};
01252 tensor H(2, h_dim, 0.0);
01253
01254 double det_of_Jacobian = 0.0;
01255
01256 tensor Jacobian;
01257
01258 double RHO;
01259 RHO= rho;
01260
01261 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
01262 {
01263 r = get_Gauss_p_c( r_integration_order, GP_c_r );
01264 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
01265 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
01266 {
01267 s = get_Gauss_p_c( s_integration_order, GP_c_s );
01268 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
01269 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
01270 {
01271 t = get_Gauss_p_c( t_integration_order, GP_c_t );
01272 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
01273
01274
01275 where =
01276 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
01277
01278 dh = dh_drst_at(r,s,t);
01279
01280 Jacobian = Jacobian_3D(dh);
01281
01282
01283
01284
01285 det_of_Jacobian = Jacobian.determinant();
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297 H = H_3D(r,s,t);
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317 weight = rw * sw * tw * RHO * det_of_Jacobian;
01318
01319 Mm = Mm + H("ib")*H("kb")*weight;
01320
01321
01322 }
01323 }
01324 }
01325
01326
01327 return Mm;
01328 }
01329
01330
01332
01333 tensor TwentySevenNodeBrick::stiffness_matrix(const tensor & K)
01334 {
01335
01336
01337 matrix Kmatrix(81,81,0.0);
01338
01339 int Ki=0;
01340 int Kj=0;
01341
01342 for ( int i=1 ; i<=27 ; i++ )
01343 {
01344 for ( int j=1 ; j<=27 ; j++ )
01345 {
01346 for ( int k=1 ; k<=3 ; k++ )
01347 {
01348 for ( int l=1 ; l<=3 ; l++ )
01349 {
01350 Ki = k+3*(i-1);
01351 Kj = l+3*(j-1);
01352
01353
01354 Kmatrix.val( Ki , Kj ) = K.cval(i,k,l,j);
01355 }
01356 }
01357 }
01358 }
01359 return Kmatrix;
01360 }
01361
01362
01363
01365 tensor TwentySevenNodeBrick::mass_matrix(const tensor & M)
01366 {
01367
01368
01369 matrix Mmatrix(81,81,0.0);
01370
01371 for ( int i=1 ; i<=81 ; i++ )
01372 {
01373 for ( int j=1 ; j<=81 ; j++ )
01374 {
01375 Mmatrix.val( i , j ) = M.cval(i,j);
01376
01377 }
01378 }
01379 return Mmatrix;
01380 }
01382
01384 tensor TwentySevenNodeBrick::Jacobian_3D(tensor dh)
01385 {
01386 tensor N_C = Nodal_Coordinates();
01387 tensor Jacobian_3D = dh("ij") * N_C("ik");
01388 return Jacobian_3D;
01389 }
01390
01391
01392 tensor TwentySevenNodeBrick::Jacobian_3Dinv(tensor dh)
01393 {
01394 tensor N_C = Nodal_Coordinates();
01395 tensor Jacobian_3Dinv = (dh("ij") * N_C("ik")).inverse();
01396 return Jacobian_3Dinv;
01397 }
01398
01399
01401 tensor TwentySevenNodeBrick::Nodal_Coordinates(void)
01402 {
01403 const int dimensions[] = {27,3};
01404 tensor N_coord(2, dimensions, 0.0);
01405
01406
01407 const Vector &nd1Crds = theNodes[0]->getCrds();
01408 const Vector &nd2Crds = theNodes[1]->getCrds();
01409 const Vector &nd3Crds = theNodes[2]->getCrds();
01410 const Vector &nd4Crds = theNodes[3]->getCrds();
01411 const Vector &nd5Crds = theNodes[4]->getCrds();
01412 const Vector &nd6Crds = theNodes[5]->getCrds();
01413 const Vector &nd7Crds = theNodes[6]->getCrds();
01414 const Vector &nd8Crds = theNodes[7]->getCrds();
01415
01416 const Vector &nd9Crds = theNodes[8]->getCrds();
01417 const Vector &nd10Crds = theNodes[9]->getCrds();
01418 const Vector &nd11Crds = theNodes[10]->getCrds();
01419 const Vector &nd12Crds = theNodes[11]->getCrds();
01420
01421 const Vector &nd13Crds = theNodes[12]->getCrds();
01422 const Vector &nd14Crds = theNodes[13]->getCrds();
01423 const Vector &nd15Crds = theNodes[14]->getCrds();
01424 const Vector &nd16Crds = theNodes[15]->getCrds();
01425
01426
01427 const Vector &nd17Crds = theNodes[16]->getCrds();
01428 const Vector &nd18Crds = theNodes[17]->getCrds();
01429 const Vector &nd19Crds = theNodes[18]->getCrds();
01430 const Vector &nd20Crds = theNodes[19]->getCrds();
01431
01432 const Vector &nd21Crds = theNodes[20]->getCrds();
01433 const Vector &nd22Crds = theNodes[21]->getCrds();
01434 const Vector &nd23Crds = theNodes[22]->getCrds();
01435 const Vector &nd24Crds = theNodes[23]->getCrds();
01436 const Vector &nd25Crds = theNodes[24]->getCrds();
01437 const Vector &nd26Crds = theNodes[25]->getCrds();
01438 const Vector &nd27Crds = theNodes[26]->getCrds();
01439
01440
01441
01442 N_coord.val(1,1)=nd1Crds(0); N_coord.val(1,2)=nd1Crds(1); N_coord.val(1,3)=nd1Crds(2);
01443 N_coord.val(2,1)=nd2Crds(0); N_coord.val(2,2)=nd2Crds(1); N_coord.val(2,3)=nd2Crds(2);
01444 N_coord.val(3,1)=nd3Crds(0); N_coord.val(3,2)=nd3Crds(1); N_coord.val(3,3)=nd3Crds(2);
01445 N_coord.val(4,1)=nd4Crds(0); N_coord.val(4,2)=nd4Crds(1); N_coord.val(4,3)=nd4Crds(2);
01446 N_coord.val(5,1)=nd5Crds(0); N_coord.val(5,2)=nd5Crds(1); N_coord.val(5,3)=nd5Crds(2);
01447 N_coord.val(6,1)=nd6Crds(0); N_coord.val(6,2)=nd6Crds(1); N_coord.val(6,3)=nd6Crds(2);
01448 N_coord.val(7,1)=nd7Crds(0); N_coord.val(7,2)=nd7Crds(1); N_coord.val(7,3)=nd7Crds(2);
01449 N_coord.val(8,1)=nd8Crds(0); N_coord.val(8,2)=nd8Crds(1); N_coord.val(8,3)=nd8Crds(2);
01450
01451 N_coord.val(9 ,1)=nd9Crds(0); N_coord.val(9 ,2)=nd9Crds(1); N_coord.val(9 ,3)=nd9Crds(2);
01452 N_coord.val(10,1)=nd10Crds(0); N_coord.val(10,2)=nd10Crds(1); N_coord.val(10,3)=nd10Crds(2);
01453 N_coord.val(11,1)=nd11Crds(0); N_coord.val(11,2)=nd11Crds(1); N_coord.val(11,3)=nd11Crds(2);
01454 N_coord.val(12,1)=nd12Crds(0); N_coord.val(12,2)=nd12Crds(1); N_coord.val(12,3)=nd12Crds(2);
01455
01456 N_coord.val(13,1)=nd13Crds(0); N_coord.val(13,2)=nd13Crds(1); N_coord.val(13,3)=nd13Crds(2);
01457 N_coord.val(14,1)=nd14Crds(0); N_coord.val(14,2)=nd14Crds(1); N_coord.val(14,3)=nd14Crds(2);
01458 N_coord.val(15,1)=nd15Crds(0); N_coord.val(15,2)=nd15Crds(1); N_coord.val(15,3)=nd15Crds(2);
01459 N_coord.val(16,1)=nd16Crds(0); N_coord.val(16,2)=nd16Crds(1); N_coord.val(16,3)=nd16Crds(2);
01460
01461 N_coord.val(17,1)=nd17Crds(0); N_coord.val(17,2)=nd17Crds(1); N_coord.val(17,3)=nd17Crds(2);
01462 N_coord.val(18,1)=nd18Crds(0); N_coord.val(18,2)=nd18Crds(1); N_coord.val(18,3)=nd18Crds(2);
01463 N_coord.val(19,1)=nd19Crds(0); N_coord.val(19,2)=nd19Crds(1); N_coord.val(19,3)=nd19Crds(2);
01464 N_coord.val(20,1)=nd20Crds(0); N_coord.val(20,2)=nd20Crds(1); N_coord.val(20,3)=nd20Crds(2);
01465
01466 N_coord.val(21,1)=nd21Crds(0); N_coord.val(21,2)=nd21Crds(1); N_coord.val(21,3)=nd21Crds(2);
01467 N_coord.val(22,1)=nd22Crds(0); N_coord.val(22,2)=nd22Crds(1); N_coord.val(22,3)=nd22Crds(2);
01468 N_coord.val(23,1)=nd23Crds(0); N_coord.val(23,2)=nd23Crds(1); N_coord.val(23,3)=nd23Crds(2);
01469 N_coord.val(24,1)=nd24Crds(0); N_coord.val(24,2)=nd24Crds(1); N_coord.val(24,3)=nd24Crds(2);
01470 N_coord.val(25,1)=nd25Crds(0); N_coord.val(25,2)=nd25Crds(1); N_coord.val(25,3)=nd25Crds(2);
01471 N_coord.val(26,1)=nd26Crds(0); N_coord.val(26,2)=nd26Crds(1); N_coord.val(26,3)=nd26Crds(2);
01472 N_coord.val(27,1)=nd27Crds(0); N_coord.val(27,2)=nd27Crds(1); N_coord.val(27,3)=nd27Crds(2);
01473
01474
01475 return N_coord;
01476
01477 }
01478
01480 tensor TwentySevenNodeBrick::incr_disp(void)
01481 {
01482 const int dimensions[] = {27,3};
01483 tensor increment_disp(2, dimensions, 0.0);
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502 const Vector &IncrDis1 = theNodes[0]->getIncrDeltaDisp();
01503 const Vector &IncrDis2 = theNodes[1]->getIncrDeltaDisp();
01504 const Vector &IncrDis3 = theNodes[2]->getIncrDeltaDisp();
01505 const Vector &IncrDis4 = theNodes[3]->getIncrDeltaDisp();
01506 const Vector &IncrDis5 = theNodes[4]->getIncrDeltaDisp();
01507 const Vector &IncrDis6 = theNodes[5]->getIncrDeltaDisp();
01508 const Vector &IncrDis7 = theNodes[6]->getIncrDeltaDisp();
01509 const Vector &IncrDis8 = theNodes[7]->getIncrDeltaDisp();
01510
01511 const Vector &IncrDis9 = theNodes[8]->getIncrDeltaDisp();
01512 const Vector &IncrDis10 = theNodes[9]->getIncrDeltaDisp();
01513 const Vector &IncrDis11 = theNodes[10]->getIncrDeltaDisp();
01514 const Vector &IncrDis12 = theNodes[11]->getIncrDeltaDisp();
01515
01516 const Vector &IncrDis13 = theNodes[12]->getIncrDeltaDisp();
01517 const Vector &IncrDis14 = theNodes[13]->getIncrDeltaDisp();
01518 const Vector &IncrDis15 = theNodes[14]->getIncrDeltaDisp();
01519 const Vector &IncrDis16 = theNodes[15]->getIncrDeltaDisp();
01520
01521 const Vector &IncrDis17 = theNodes[16]->getIncrDeltaDisp();
01522 const Vector &IncrDis18 = theNodes[17]->getIncrDeltaDisp();
01523 const Vector &IncrDis19 = theNodes[18]->getIncrDeltaDisp();
01524 const Vector &IncrDis20 = theNodes[19]->getIncrDeltaDisp();
01525
01526 const Vector &IncrDis21 = theNodes[20]->getIncrDeltaDisp();
01527 const Vector &IncrDis22 = theNodes[21]->getIncrDeltaDisp();
01528 const Vector &IncrDis23 = theNodes[22]->getIncrDeltaDisp();
01529 const Vector &IncrDis24 = theNodes[23]->getIncrDeltaDisp();
01530 const Vector &IncrDis25 = theNodes[24]->getIncrDeltaDisp();
01531 const Vector &IncrDis26 = theNodes[25]->getIncrDeltaDisp();
01532 const Vector &IncrDis27 = theNodes[26]->getIncrDeltaDisp();
01533
01534
01535
01536 increment_disp.val(1,1)=IncrDis1(0); increment_disp.val(1,2)=IncrDis1(1);increment_disp.val(1,3)=IncrDis1(2);
01537 increment_disp.val(2,1)=IncrDis2(0); increment_disp.val(2,2)=IncrDis2(1);increment_disp.val(2,3)=IncrDis2(2);
01538 increment_disp.val(3,1)=IncrDis3(0); increment_disp.val(3,2)=IncrDis3(1);increment_disp.val(3,3)=IncrDis3(2);
01539 increment_disp.val(4,1)=IncrDis4(0); increment_disp.val(4,2)=IncrDis4(1);increment_disp.val(4,3)=IncrDis4(2);
01540 increment_disp.val(5,1)=IncrDis5(0); increment_disp.val(5,2)=IncrDis5(1);increment_disp.val(5,3)=IncrDis5(2);
01541 increment_disp.val(6,1)=IncrDis6(0); increment_disp.val(6,2)=IncrDis6(1);increment_disp.val(6,3)=IncrDis6(2);
01542 increment_disp.val(7,1)=IncrDis7(0); increment_disp.val(7,2)=IncrDis7(1);increment_disp.val(7,3)=IncrDis7(2);
01543 increment_disp.val(8,1)=IncrDis8(0); increment_disp.val(8,2)=IncrDis8(1);increment_disp.val(8,3)=IncrDis8(2);
01544
01545 increment_disp.val(9 ,1)=IncrDis9(0); increment_disp.val(9 ,2)=IncrDis9(1); increment_disp.val(9 ,3)=IncrDis9(2);
01546 increment_disp.val(10,1)=IncrDis10(0); increment_disp.val(10,2)=IncrDis10(1);increment_disp.val(10,3)=IncrDis10(2);
01547 increment_disp.val(11,1)=IncrDis11(0); increment_disp.val(11,2)=IncrDis11(1);increment_disp.val(11,3)=IncrDis11(2);
01548 increment_disp.val(12,1)=IncrDis12(0); increment_disp.val(12,2)=IncrDis12(1);increment_disp.val(12,3)=IncrDis12(2);
01549
01550 increment_disp.val(13,1)=IncrDis13(0); increment_disp.val(13,2)=IncrDis13(1);increment_disp.val(13,3)=IncrDis13(2);
01551 increment_disp.val(14,1)=IncrDis14(0); increment_disp.val(14,2)=IncrDis14(1);increment_disp.val(14,3)=IncrDis14(2);
01552 increment_disp.val(15,1)=IncrDis15(0); increment_disp.val(15,2)=IncrDis15(1);increment_disp.val(15,3)=IncrDis15(2);
01553 increment_disp.val(16,1)=IncrDis16(0); increment_disp.val(16,2)=IncrDis16(1);increment_disp.val(16,3)=IncrDis16(2);
01554
01555 increment_disp.val(17,1)=IncrDis17(0); increment_disp.val(17,2)=IncrDis17(1);increment_disp.val(17,3)=IncrDis17(2);
01556 increment_disp.val(18,1)=IncrDis18(0); increment_disp.val(18,2)=IncrDis18(1);increment_disp.val(18,3)=IncrDis18(2);
01557 increment_disp.val(19,1)=IncrDis19(0); increment_disp.val(19,2)=IncrDis19(1);increment_disp.val(19,3)=IncrDis19(2);
01558 increment_disp.val(20,1)=IncrDis20(0); increment_disp.val(20,2)=IncrDis20(1);increment_disp.val(20,3)=IncrDis20(2);
01559
01560 increment_disp.val(21,1)=IncrDis21(0); increment_disp.val(21,2)=IncrDis21(1);increment_disp.val(21,3)=IncrDis21(2);
01561 increment_disp.val(22,1)=IncrDis22(0); increment_disp.val(22,2)=IncrDis22(1);increment_disp.val(22,3)=IncrDis22(2);
01562 increment_disp.val(23,1)=IncrDis23(0); increment_disp.val(23,2)=IncrDis23(1);increment_disp.val(23,3)=IncrDis23(2);
01563 increment_disp.val(24,1)=IncrDis24(0); increment_disp.val(24,2)=IncrDis24(1);increment_disp.val(24,3)=IncrDis24(2);
01564 increment_disp.val(25,1)=IncrDis25(0); increment_disp.val(25,2)=IncrDis25(1);increment_disp.val(25,3)=IncrDis25(2);
01565 increment_disp.val(26,1)=IncrDis26(0); increment_disp.val(26,2)=IncrDis26(1);increment_disp.val(26,3)=IncrDis26(2);
01566 increment_disp.val(27,1)=IncrDis27(0); increment_disp.val(27,2)=IncrDis27(1);increment_disp.val(27,3)=IncrDis27(2);
01567
01568
01569
01570 return increment_disp;
01571 }
01572
01574 tensor TwentySevenNodeBrick::total_disp(void)
01575 {
01576 const int dimensions[] = {27,3};
01577 tensor total_disp(2, dimensions, 0.0);
01578
01579
01580 const Vector &TotDis1 = theNodes[0]->getTrialDisp();
01581 opserr<<"\ntot node " << theNodes[0]->getTag() <<" x "<< TotDis1(0) <<" y "<< TotDis1(1) << " z "<< TotDis1(2) << endln;
01582 const Vector &TotDis2 = theNodes[1]->getTrialDisp();
01583 opserr << "tot node " << theNodes[1]->getTag() << " x " << TotDis2(0) <<" y "<< TotDis2(1) << " z "<< TotDis2(2) << endln;
01584 const Vector &TotDis3 = theNodes[2]->getTrialDisp();
01585 opserr << "tot node " << theNodes[2]->getTag() << " x " << TotDis3(0) <<" y "<< TotDis3(1) << " z "<< TotDis3(2) << endln;
01586 const Vector &TotDis4 = theNodes[3]->getTrialDisp();
01587 opserr << "tot node " << theNodes[3]->getTag() << " x " << TotDis4(0) <<" y "<< TotDis4(1) << " z "<< TotDis4(2) << endln;
01588 const Vector &TotDis5 = theNodes[4]->getTrialDisp();
01589 opserr << "tot node " << theNodes[4]->getTag() << " x " << TotDis5(0) <<" y "<< TotDis5(1) << " z "<< TotDis5(2) << endln;
01590 const Vector &TotDis6 = theNodes[5]->getTrialDisp();
01591 opserr << "tot node " << theNodes[5]->getTag() << " x " << TotDis6(0) <<" y "<< TotDis6(1) << " z "<< TotDis6(2) << endln;
01592 const Vector &TotDis7 = theNodes[6]->getTrialDisp();
01593 opserr << "tot node " << theNodes[6]->getTag() << " x " << TotDis7(0) <<" y "<< TotDis7(1) << " z "<< TotDis7(2) << endln;
01594 const Vector &TotDis8 = theNodes[7]->getTrialDisp();
01595 opserr << "tot node " << theNodes[7]->getTag() << " x " << TotDis8(0) <<" y "<< TotDis8(1) << " z "<< TotDis8(2) << endln;
01596
01597 const Vector &TotDis9 = theNodes[8]->getTrialDisp();
01598 opserr << "tot node " << theNodes[8]->getTag() << " x " << TotDis9(0) <<" y "<< TotDis9(1) << " z "<< TotDis9(2) << endln;
01599 const Vector &TotDis10 = theNodes[9]->getTrialDisp();
01600 opserr << "tot node " << theNodes[9]->getTag() << " x " << TotDis10(0) <<" y "<< TotDis10(1) << " z "<< TotDis10(2) << endln;
01601 const Vector &TotDis11 = theNodes[10]->getTrialDisp();
01602 opserr << "tot node " << theNodes[10]->getTag() << " x " << TotDis11(0) <<" y "<< TotDis11(1) << " z "<< TotDis11(2) << endln;
01603 const Vector &TotDis12 = theNodes[11]->getTrialDisp();
01604 opserr << "tot node " << theNodes[11]->getTag() << " x " << TotDis12(0) <<" y "<< TotDis12(1) << " z "<< TotDis12(2) << endln;
01605
01606 const Vector &TotDis13 = theNodes[12]->getTrialDisp();
01607 opserr << "tot node " << theNodes[12]->getTag() << " x " << TotDis13(0) <<" y "<< TotDis13(1) << " z "<< TotDis13(2) << endln;
01608 const Vector &TotDis14 = theNodes[13]->getTrialDisp();
01609 opserr << "tot node " << theNodes[13]->getTag() << " x " << TotDis14(0) <<" y "<< TotDis14(1) << " z "<< TotDis14(2) << endln;
01610 const Vector &TotDis15 = theNodes[14]->getTrialDisp();
01611 opserr << "tot node " << theNodes[14]->getTag() << " x " << TotDis15(0) <<" y "<< TotDis15(1) << " z "<< TotDis15(2) << endln;
01612 const Vector &TotDis16 = theNodes[15]->getTrialDisp();
01613 opserr << "tot node " << theNodes[15]->getTag() << " x " << TotDis16(0) <<" y "<< TotDis16(1) << " z "<< TotDis16(2) << endln;
01614
01615 const Vector &TotDis17 = theNodes[16]->getTrialDisp();
01616 opserr << "tot node " << theNodes[16]->getTag() << " x " << TotDis17(0) <<" y "<< TotDis17(1) << " z "<< TotDis17(2) << endln;
01617 const Vector &TotDis18 = theNodes[17]->getTrialDisp();
01618 opserr << "tot node " << theNodes[17]->getTag() << " x " << TotDis18(0) <<" y "<< TotDis18(1) << " z "<< TotDis18(2) << endln;
01619 const Vector &TotDis19 = theNodes[18]->getTrialDisp();
01620 opserr << "tot node " << theNodes[18]->getTag() << " x " << TotDis19(0) <<" y "<< TotDis19(1) << " z "<< TotDis19(2) << endln;
01621 const Vector &TotDis20 = theNodes[19]->getTrialDisp();
01622 opserr << "tot node " << theNodes[19]->getTag() << " x " << TotDis20(0) <<" y "<< TotDis20(1) << " z "<< TotDis20(2) << endln;
01623 const Vector &TotDis21 = theNodes[20]->getTrialDisp();
01624 opserr << "tot node " << theNodes[20]->getTag() << " x " << TotDis21(0) <<" y "<< TotDis21(1) << " z "<< TotDis21(2) << endln;
01625 const Vector &TotDis22 = theNodes[21]->getTrialDisp();
01626 opserr << "tot node " << theNodes[21]->getTag() << " x " << TotDis22(0) <<" y "<< TotDis22(1) << " z "<< TotDis22(2) << endln;
01627 const Vector &TotDis23 = theNodes[22]->getTrialDisp();
01628 opserr << "tot node " << theNodes[22]->getTag() << " x " << TotDis23(0) <<" y "<< TotDis23(1) << " z "<< TotDis23(2) << endln;
01629 const Vector &TotDis24 = theNodes[23]->getTrialDisp();
01630 opserr << "tot node " << theNodes[23]->getTag() << " x " << TotDis24(0) <<" y "<< TotDis24(1) << " z "<< TotDis24(2) << endln;
01631 const Vector &TotDis25 = theNodes[24]->getTrialDisp();
01632 opserr << "tot node " << theNodes[24]->getTag() << " x " << TotDis25(0) <<" y "<< TotDis25(1) << " z "<< TotDis25(2) << endln;
01633 const Vector &TotDis26 = theNodes[25]->getTrialDisp();
01634 opserr << "tot node " << theNodes[25]->getTag() << " x " << TotDis26(0) <<" y "<< TotDis26(1) << " z "<< TotDis26(2) << endln;
01635 const Vector &TotDis27 = theNodes[26]->getTrialDisp();
01636 opserr << "tot node " << theNodes[26]->getTag() << " x " << TotDis27(0) <<" y "<< TotDis27(1) << " z "<< TotDis27(2) << endln;
01637
01638
01639
01640
01641 total_disp.val(1,1)=TotDis1(0); total_disp.val(1,2)=TotDis1(1);total_disp.val(1,3)=TotDis1(2);
01642 total_disp.val(2,1)=TotDis2(0); total_disp.val(2,2)=TotDis2(1);total_disp.val(2,3)=TotDis2(2);
01643 total_disp.val(3,1)=TotDis3(0); total_disp.val(3,2)=TotDis3(1);total_disp.val(3,3)=TotDis3(2);
01644 total_disp.val(4,1)=TotDis4(0); total_disp.val(4,2)=TotDis4(1);total_disp.val(4,3)=TotDis4(2);
01645 total_disp.val(5,1)=TotDis5(0); total_disp.val(5,2)=TotDis5(1);total_disp.val(5,3)=TotDis5(2);
01646 total_disp.val(6,1)=TotDis6(0); total_disp.val(6,2)=TotDis6(1);total_disp.val(6,3)=TotDis6(2);
01647 total_disp.val(7,1)=TotDis7(0); total_disp.val(7,2)=TotDis7(1);total_disp.val(7,3)=TotDis7(2);
01648 total_disp.val(8,1)=TotDis8(0); total_disp.val(8,2)=TotDis8(1);total_disp.val(8,3)=TotDis8(2);
01649
01650 total_disp.val(9,1)=TotDis9(0); total_disp.val(9,2)=TotDis9(1);total_disp.val(9,3)=TotDis9(2);
01651 total_disp.val(10,1)=TotDis10(0); total_disp.val(10,2)=TotDis10(1);total_disp.val(10,3)=TotDis10(2);
01652 total_disp.val(11,1)=TotDis11(0); total_disp.val(11,2)=TotDis11(1);total_disp.val(11,3)=TotDis11(2);
01653 total_disp.val(12,1)=TotDis12(0); total_disp.val(12,2)=TotDis12(1);total_disp.val(12,3)=TotDis12(2);
01654
01655 total_disp.val(13,1)=TotDis13(0); total_disp.val(13,2)=TotDis13(1);total_disp.val(13,3)=TotDis13(2);
01656 total_disp.val(14,1)=TotDis14(0); total_disp.val(14,2)=TotDis14(1);total_disp.val(14,3)=TotDis14(2);
01657 total_disp.val(15,1)=TotDis15(0); total_disp.val(15,2)=TotDis15(1);total_disp.val(15,3)=TotDis15(2);
01658 total_disp.val(16,1)=TotDis16(0); total_disp.val(16,2)=TotDis16(1);total_disp.val(16,3)=TotDis16(2);
01659
01660 total_disp.val(17,1)=TotDis17(0); total_disp.val(17,2)=TotDis17(1);total_disp.val(17,3)=TotDis17(2);
01661 total_disp.val(18,1)=TotDis18(0); total_disp.val(18,2)=TotDis18(1);total_disp.val(18,3)=TotDis18(2);
01662 total_disp.val(19,1)=TotDis19(0); total_disp.val(19,2)=TotDis19(1);total_disp.val(19,3)=TotDis19(2);
01663 total_disp.val(20,1)=TotDis20(0); total_disp.val(20,2)=TotDis20(1);total_disp.val(20,3)=TotDis20(2);
01664
01665 total_disp.val(21,1)=TotDis21(0); total_disp.val(21,2)=TotDis21(1);total_disp.val(21,3)=TotDis21(2);
01666 total_disp.val(22,1)=TotDis22(0); total_disp.val(22,2)=TotDis22(1);total_disp.val(22,3)=TotDis22(2);
01667 total_disp.val(23,1)=TotDis23(0); total_disp.val(23,2)=TotDis23(1);total_disp.val(23,3)=TotDis23(2);
01668 total_disp.val(24,1)=TotDis24(0); total_disp.val(24,2)=TotDis24(1);total_disp.val(24,3)=TotDis24(2);
01669 total_disp.val(25,1)=TotDis25(0); total_disp.val(25,2)=TotDis25(1);total_disp.val(25,3)=TotDis25(2);
01670 total_disp.val(26,1)=TotDis26(0); total_disp.val(26,2)=TotDis26(1);total_disp.val(26,3)=TotDis26(2);
01671 total_disp.val(27,1)=TotDis27(0); total_disp.val(27,2)=TotDis27(1);total_disp.val(27,3)=TotDis27(2);
01672
01673
01674 return total_disp;
01675 }
01676
01677
01679 tensor TwentySevenNodeBrick::total_disp(FILE *fp, double * u)
01680 {
01681 const int dimensions[] = {27,3};
01682 tensor total_disp(2, dimensions, 0.0);
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705 const Vector &TotDis1 = theNodes[0]->getTrialDisp();
01706 const Vector &TotDis2 = theNodes[1]->getTrialDisp();
01707 const Vector &TotDis3 = theNodes[2]->getTrialDisp();
01708 const Vector &TotDis4 = theNodes[3]->getTrialDisp();
01709 const Vector &TotDis5 = theNodes[4]->getTrialDisp();
01710 const Vector &TotDis6 = theNodes[5]->getTrialDisp();
01711 const Vector &TotDis7 = theNodes[6]->getTrialDisp();
01712 const Vector &TotDis8 = theNodes[7]->getTrialDisp();
01713
01714 total_disp.val(1,1)=TotDis1(0); total_disp.val(1,2)=TotDis1(1);total_disp.val(1,3)=TotDis1(2);
01715 total_disp.val(2,1)=TotDis2(0); total_disp.val(2,2)=TotDis2(1);total_disp.val(2,3)=TotDis2(2);
01716 total_disp.val(3,1)=TotDis3(0); total_disp.val(3,2)=TotDis3(1);total_disp.val(3,3)=TotDis3(2);
01717 total_disp.val(4,1)=TotDis4(0); total_disp.val(4,2)=TotDis4(1);total_disp.val(4,3)=TotDis4(2);
01718 total_disp.val(5,1)=TotDis5(0); total_disp.val(5,2)=TotDis5(1);total_disp.val(5,3)=TotDis5(2);
01719 total_disp.val(6,1)=TotDis6(0); total_disp.val(6,2)=TotDis6(1);total_disp.val(6,3)=TotDis6(2);
01720 total_disp.val(7,1)=TotDis7(0); total_disp.val(7,2)=TotDis7(1);total_disp.val(7,3)=TotDis7(2);
01721 total_disp.val(8,1)=TotDis8(0); total_disp.val(8,2)=TotDis8(1);total_disp.val(8,3)=TotDis8(2);
01722
01723 return total_disp;
01724 }
01725
01726
01728 int TwentySevenNodeBrick::get_global_number_of_node(int local_node_number)
01729 {
01730
01731 return connectedExternalNodes(local_node_number);
01732 }
01733
01735 int TwentySevenNodeBrick::get_Brick_Number(void)
01736 {
01737
01738 return this->getTag();
01739 }
01740
01742
01743
01744
01745
01746
01747
01748
01750
01751
01756
01757
01759
01760
01761
01762
01763
01764
01765
01766
01772
01773
01774
01775
01777
01778 tensor TwentySevenNodeBrick::nodal_forces(void)
01779 {
01780 int force_dim[] = {27,3};
01781
01782 tensor nodal_forces(2,force_dim,0.0);
01783
01784 double r = 0.0;
01785 double rw = 0.0;
01786 double s = 0.0;
01787 double sw = 0.0;
01788 double t = 0.0;
01789 double tw = 0.0;
01790
01791 short where = 0;
01792 double weight = 0.0;
01793
01794 int dh_dim[] = {27,3};
01795
01796 tensor dh(2, dh_dim, 0.0);
01797
01798 stresstensor stress_at_GP(0.0);
01799
01800 double det_of_Jacobian = 0.0;
01801
01802 straintensor incremental_strain;
01803
01804 static int disp_dim[] = {27,3};
01805 tensor incremental_displacements(2,disp_dim,0.0);
01806
01807 incremental_displacements = incr_disp();
01808
01809 tensor Jacobian;
01810 tensor JacobianINV;
01811 tensor dhGlobal;
01812
01813 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
01814 {
01815 r = get_Gauss_p_c( r_integration_order, GP_c_r );
01816 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
01817
01818 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
01819 {
01820 s = get_Gauss_p_c( s_integration_order, GP_c_s );
01821 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
01822
01823 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
01824 {
01825 t = get_Gauss_p_c( t_integration_order, GP_c_t );
01826 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
01827
01828
01829
01830 where =
01831 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
01832
01833
01834 dh = dh_drst_at(r,s,t);
01835
01836
01837 Jacobian = Jacobian_3D(dh);
01838
01839
01840
01841 JacobianINV = Jacobian_3Dinv(dh);
01842
01843
01844
01845 det_of_Jacobian = Jacobian.determinant();
01846
01847
01848
01849
01850 dhGlobal = dh("ij") * JacobianINV("kj");
01851
01852
01853 weight = rw * sw * tw * det_of_Jacobian;
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01866
01867
01869
01870
01871
01872
01874
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928 stress_at_GP = matpoint[where]->getStressTensor();
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962 nodal_forces =
01963 nodal_forces + dhGlobal("ib")*stress_at_GP("ab")*weight;
01964
01965
01966 }
01967 }
01968 }
01969
01970
01971
01972 return nodal_forces;
01973
01974 }
01975
01977
01978 tensor TwentySevenNodeBrick::iterative_nodal_forces(void)
01979 {
01980 int force_dim[] = {27,3};
01981
01982 tensor nodal_forces(2,force_dim,0.0);
01983
01984 double r = 0.0;
01985 double rw = 0.0;
01986 double s = 0.0;
01987 double sw = 0.0;
01988 double t = 0.0;
01989 double tw = 0.0;
01990
01991 short where = 0;
01992 double weight = 0.0;
01993
01994 int dh_dim[] = {27,3};
01995
01996 tensor dh(2, dh_dim, 0.0);
01997
01998 stresstensor stress_at_GP(0.0);
01999
02000 double det_of_Jacobian = 0.0;
02001
02002 tensor Jacobian;
02003 tensor JacobianINV;
02004 tensor dhGlobal;
02005
02006 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02007 {
02008 r = get_Gauss_p_c( r_integration_order, GP_c_r );
02009 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
02010
02011 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02012 {
02013 s = get_Gauss_p_c( s_integration_order, GP_c_s );
02014 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
02015
02016 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02017 {
02018 t = get_Gauss_p_c( t_integration_order, GP_c_t );
02019 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
02020
02021
02022
02023 where =
02024 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
02025
02026
02027
02028
02029
02030
02031 dh = dh_drst_at(r,s,t);
02032
02033
02034 Jacobian = Jacobian_3D(dh);
02035
02036
02037
02038 JacobianINV = Jacobian_3Dinv(dh);
02039
02040
02041
02042 det_of_Jacobian = Jacobian.determinant();
02043
02044
02045
02046
02047 dhGlobal = dh("ij") * JacobianINV("kj");
02048
02049
02050 weight = rw * sw * tw * det_of_Jacobian;
02051
02052
02053
02054
02055
02056 stress_at_GP = matpoint[where]->getStressTensor();
02057 stress_at_GP.reportshortpqtheta("\n iterative_stress at GAUSS point in iterative_nodal_force\n");
02058
02059
02060 nodal_forces =
02061 nodal_forces + dhGlobal("ib")*stress_at_GP("ab")*weight;
02062
02063
02064 }
02065 }
02066 }
02067
02068
02069 return nodal_forces;
02070
02071 }
02072
02074
02075 tensor TwentySevenNodeBrick::nodal_forces_from_stress(stresstensor & stress)
02076 {
02077 int force_dim[] = {27,3};
02078
02079 tensor nodal_forces(2,force_dim,0.0);
02080
02081 double r = 0.0;
02082 double rw = 0.0;
02083 double s = 0.0;
02084 double sw = 0.0;
02085 double t = 0.0;
02086 double tw = 0.0;
02087
02088 double weight = 0.0;
02089
02090 int dh_dim[] = {27,3};
02091
02092 tensor dh(2, dh_dim, 0.0);
02093
02094 double det_of_Jacobian = 0.0;
02095
02096 tensor Jacobian;
02097 tensor JacobianINV;
02098 tensor dhGlobal;
02099
02100 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02101 {
02102 r = get_Gauss_p_c( r_integration_order, GP_c_r );
02103 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
02104
02105 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02106 {
02107 s = get_Gauss_p_c( s_integration_order, GP_c_s );
02108 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
02109
02110 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02111 {
02112 t = get_Gauss_p_c( t_integration_order, GP_c_t );
02113 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125 dh = dh_drst_at(r,s,t);
02126
02127
02128 Jacobian = Jacobian_3D(dh);
02129
02130
02131
02132 JacobianINV = Jacobian_3Dinv(dh);
02133
02134
02135
02136 det_of_Jacobian = Jacobian.determinant();
02137
02138
02139
02140
02141 dhGlobal = dh("ij") * JacobianINV("kj");
02142
02143
02144 weight = rw * sw * tw * det_of_Jacobian;
02145
02146
02147
02148
02149
02150
02151 nodal_forces =
02152 nodal_forces + dhGlobal("ib")*stress("ab")*weight;
02153
02154
02155 }
02156 }
02157 }
02158
02159 return nodal_forces;
02160
02161 }
02162
02163
02165
02166
02167 tensor TwentySevenNodeBrick::linearized_nodal_forces(void)
02168 {
02169 int force_dim[] = {27,3};
02170
02171 tensor linearized_nodal_forces(2,force_dim,0.0);
02172
02173 double r = 0.0;
02174 double rw = 0.0;
02175 double s = 0.0;
02176 double sw = 0.0;
02177 double t = 0.0;
02178 double tw = 0.0;
02179
02180 short where = 0;
02181 double weight = 0.0;
02182
02183 int dh_dim[] = {27,3};
02184
02185 tensor dh(2, dh_dim, 0.0);
02186
02187 tensor Constitutive( 4, def_dim_4, 0.0);
02188
02189 double det_of_Jacobian = 0.0;
02190
02191 static int disp_dim[] = {27,3};
02192
02193 tensor incremental_displacements(2,disp_dim,0.0);
02194
02195 straintensor incremental_strain;
02196
02197 tensor Jacobian;
02198 tensor JacobianINV;
02199 tensor dhGlobal;
02200
02201 stresstensor final_linearized_stress;
02202
02203
02204 incremental_displacements = incr_disp();
02205
02206
02207 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02208 {
02209 r = get_Gauss_p_c( r_integration_order, GP_c_r );
02210 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
02211
02212 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02213 {
02214 s = get_Gauss_p_c( s_integration_order, GP_c_s );
02215 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
02216
02217 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02218 {
02219 t = get_Gauss_p_c( t_integration_order, GP_c_t );
02220 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
02221
02222
02223
02224 where =
02225 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
02226
02227
02228 dh = dh_drst_at(r,s,t);
02229
02230
02231 Jacobian = Jacobian_3D(dh);
02232
02233
02234
02235 JacobianINV = Jacobian_3Dinv(dh);
02236
02237
02238
02239 det_of_Jacobian = Jacobian.determinant();
02240
02241
02242
02243
02244 dhGlobal = dh("ij") * JacobianINV("kj");
02245
02246
02247 weight = rw * sw * tw * det_of_Jacobian;
02248
02249
02250
02251
02252
02253
02254
02255
02256 incremental_strain =
02257 (dhGlobal("ib")*incremental_displacements("ia")).symmetrize11();
02258 incremental_strain.null_indices();
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268 if ( ! (matpoint[where]->matmodel)->setTrialStrainIncr( incremental_strain) )
02269 opserr << "TwentySevenNodeBrick::linearized_nodal_forces (tag: " << this->getTag() << "), not converged\n";
02270
02271 Constitutive = (matpoint[where]->matmodel)->getTangentTensor();
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286 final_linearized_stress =
02287 Constitutive("ijkl") * incremental_strain("kl");
02288
02289
02290 linearized_nodal_forces = linearized_nodal_forces +
02291 dhGlobal("ib")*final_linearized_stress("ab")*weight;
02292
02293
02294 }
02295 }
02296 }
02297
02298
02299 return linearized_nodal_forces;
02300
02301 }
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02409
02410
02411
02412
02413
02414
02415
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427 void TwentySevenNodeBrick::report(char * msg)
02428 {
02429 if ( msg ) ::printf("** %s",msg);
02430 ::printf("\n Element Number = %d\n", this->getTag() );
02431 ::printf("\n Number of nodes in a TwentySevenNodebrick = %d\n",
02432 nodes_in_brick);
02433 ::printf("\n Determinant of Jacobian (! ==0 before comp.) = %f\n",
02434 determinant_of_Jacobian);
02435
02436 ::printf("Node numbers \n");
02437 ::printf(".....1.....2.....3.....4.....5.....6.....7.....8.....9.....0.....1.....2\n");
02438 for ( int i=0 ; i<27 ; i++ )
02439
02440 ::printf("%6d",connectedExternalNodes(i));
02441 ::printf("\n");
02442
02443
02444 ::printf("\n\n");
02445
02446
02447
02448
02449
02450
02451
02452 int total_number_of_Gauss_points = r_integration_order*
02453 s_integration_order*
02454 t_integration_order;
02455 if ( total_number_of_Gauss_points != 0 )
02456 {
02457
02458
02459
02460
02461
02462
02463 theNodes[0]->Print(opserr);
02464 theNodes[1]->Print(opserr);
02465 theNodes[2]->Print(opserr);
02466 theNodes[3]->Print(opserr);
02467 theNodes[4]->Print(opserr);
02468 theNodes[5]->Print(opserr);
02469 theNodes[6]->Print(opserr);
02470 theNodes[7]->Print(opserr);
02471 theNodes[8]->Print(opserr);
02472 theNodes[9]->Print(opserr);
02473 theNodes[10]->Print(opserr);
02474 theNodes[11]->Print(opserr);
02475 theNodes[12]->Print(opserr);
02476 theNodes[13]->Print(opserr);
02477 theNodes[14]->Print(opserr);
02478 theNodes[15]->Print(opserr);
02479 theNodes[16]->Print(opserr);
02480 theNodes[17]->Print(opserr);
02481 theNodes[18]->Print(opserr);
02482 theNodes[19]->Print(opserr);
02483 theNodes[20]->Print(opserr);
02484 theNodes[21]->Print(opserr);
02485 theNodes[22]->Print(opserr);
02486 theNodes[23]->Print(opserr);
02487 theNodes[24]->Print(opserr);
02488 theNodes[25]->Print(opserr);
02489 theNodes[26]->Print(opserr);
02490
02491
02492
02493
02494
02495
02496
02497 }
02498
02499 ::printf("\n\nGauss-Legendre integration order\n");
02500 ::printf("Integration order in r direction = %d\n",r_integration_order);
02501 ::printf("Integration order in s direction = %d\n",s_integration_order);
02502 ::printf("Integration order in t direction = %d\n\n",t_integration_order);
02503
02504
02505
02506 for( int GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02507 {
02508 for( int GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02509 {
02510 for( int GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02511 {
02512
02513
02514 short where =
02515 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
02516
02517 ::printf("\n\n----------------**************** where = %d \n", where);
02518 ::printf(" GP_c_r = %d, GP_c_s = %d, GP_c_t = %d\n",
02519 GP_c_r,GP_c_s,GP_c_t);
02520 matpoint[where]->report("Material Point\n");
02521
02522
02523
02524 }
02525 }
02526 }
02527
02528 }
02529
02530
02531
02532 void TwentySevenNodeBrick::reportshort(char * msg)
02533 {
02534 if ( msg ) ::printf("** %s",msg);
02535 ::printf("\n Element Number = %d\n", this->getTag() );
02536 ::printf("\n Number of nodes in a TwentySevenNodeBrick = %d\n",
02537 nodes_in_brick);
02538 ::printf("\n Determinant of Jacobian (! ==0 before comp.) = %f\n",
02539 determinant_of_Jacobian);
02540
02541 ::printf("Node numbers \n");
02542 ::printf(".....1.....2.....3.....4.....5.....6.....7.....8.....9.....0.....1.....2\n");
02543 for ( int i=0 ; i<nodes_in_brick ; i++ )
02544
02545 ::printf( "%6d",connectedExternalNodes(i) );
02546
02547 ::printf("\n");
02548
02549
02550 ::printf("\n\n");
02551
02552
02553
02554
02555 ::printf("\n\n");
02556
02557 }
02558
02559
02560
02561
02562
02563 void TwentySevenNodeBrick::reportPAK(char * msg)
02564 {
02565 if ( msg ) ::printf("%s",msg);
02566 ::printf("%10d ", this->getTag());
02567 for ( int i=0 ; i<nodes_in_brick ; i++ )
02568 ::printf( "%6d",connectedExternalNodes(i) );
02569
02570
02571 printf("\n");
02572 }
02573
02574
02575
02576 void TwentySevenNodeBrick::reportpqtheta(int GP_numb)
02577 {
02578 short where = GP_numb-1;
02579 matpoint[where]->reportpqtheta("");
02580 }
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597 void TwentySevenNodeBrick::computeGaussPoint(void)
02598 {
02599
02600
02601
02602
02603 int count;
02604 count = FixedOrder*FixedOrder*FixedOrder;
02605
02606 GaussCoord(0) = count;
02607
02608 double r = 0.0;
02609 double s = 0.0;
02610 double t = 0.0;
02611
02612 short where = 0;
02613
02614
02615 static const int dim[] = {3, 27};
02616 static const int dim27[] = {3, count};
02617 tensor NodalCoord(2, dim, 0.0);
02618 tensor matpointCoord(2, dim27, 0.0);
02619 int h_dim[] = {81,3};
02620 tensor H(2, h_dim, 0.0);
02621
02622
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641
02642
02643
02644 const Vector &nd1Crds = theNodes[0]->getCrds();
02645 const Vector &nd2Crds = theNodes[1]->getCrds();
02646 const Vector &nd3Crds = theNodes[2]->getCrds();
02647 const Vector &nd4Crds = theNodes[3]->getCrds();
02648
02649 const Vector &nd5Crds = theNodes[4]->getCrds();
02650 const Vector &nd6Crds = theNodes[5]->getCrds();
02651 const Vector &nd7Crds = theNodes[6]->getCrds();
02652 const Vector &nd8Crds = theNodes[7]->getCrds();
02653
02654 const Vector &nd9Crds = theNodes[8]->getCrds();
02655 const Vector &nd10Crds = theNodes[9]->getCrds();
02656 const Vector &nd11Crds = theNodes[10]->getCrds();
02657 const Vector &nd12Crds = theNodes[11]->getCrds();
02658
02659 const Vector &nd13Crds = theNodes[12]->getCrds();
02660 const Vector &nd14Crds = theNodes[13]->getCrds();
02661 const Vector &nd15Crds = theNodes[14]->getCrds();
02662 const Vector &nd16Crds = theNodes[15]->getCrds();
02663
02664 const Vector &nd17Crds = theNodes[16]->getCrds();
02665 const Vector &nd18Crds = theNodes[17]->getCrds();
02666 const Vector &nd19Crds = theNodes[18]->getCrds();
02667 const Vector &nd20Crds = theNodes[19]->getCrds();
02668
02669 const Vector &nd21Crds = theNodes[20]->getCrds();
02670 const Vector &nd22Crds = theNodes[21]->getCrds();
02671 const Vector &nd23Crds = theNodes[22]->getCrds();
02672 const Vector &nd24Crds = theNodes[23]->getCrds();
02673 const Vector &nd25Crds = theNodes[24]->getCrds();
02674 const Vector &nd26Crds = theNodes[25]->getCrds();
02675 const Vector &nd27Crds = theNodes[26]->getCrds();
02676
02677
02678
02679 NodalCoord.val(1, 1)=nd1Crds( 0); NodalCoord.val(2, 1)=nd1Crds( 1); NodalCoord.val(3, 1)=nd1Crds( 2);
02680 NodalCoord.val(1, 2)=nd2Crds( 0); NodalCoord.val(2, 2)=nd2Crds( 1); NodalCoord.val(3, 2)=nd2Crds( 2);
02681 NodalCoord.val(1, 3)=nd3Crds( 0); NodalCoord.val(2, 3)=nd3Crds( 1); NodalCoord.val(3, 3)=nd3Crds( 2);
02682 NodalCoord.val(1, 4)=nd4Crds( 0); NodalCoord.val(2, 4)=nd4Crds( 1); NodalCoord.val(3, 4)=nd4Crds( 2);
02683 NodalCoord.val(1, 5)=nd5Crds( 0); NodalCoord.val(2, 5)=nd5Crds( 1); NodalCoord.val(3, 5)=nd5Crds( 2);
02684 NodalCoord.val(1, 6)=nd6Crds( 0); NodalCoord.val(2, 6)=nd6Crds( 1); NodalCoord.val(3, 6)=nd6Crds( 2);
02685 NodalCoord.val(1, 7)=nd7Crds( 0); NodalCoord.val(2, 7)=nd7Crds( 1); NodalCoord.val(3, 7)=nd7Crds( 2);
02686 NodalCoord.val(1, 8)=nd8Crds( 0); NodalCoord.val(2, 8)=nd8Crds( 1); NodalCoord.val(3, 8)=nd8Crds( 2);
02687 NodalCoord.val(1, 9)=nd9Crds( 0); NodalCoord.val(2, 9)=nd9Crds( 1); NodalCoord.val(3, 9)=nd9Crds( 2);
02688 NodalCoord.val(1,10)=nd10Crds(0); NodalCoord.val(2,10)=nd10Crds(1); NodalCoord.val(3,10)=nd10Crds(2);
02689 NodalCoord.val(1,11)=nd11Crds(0); NodalCoord.val(2,11)=nd11Crds(1); NodalCoord.val(3,11)=nd11Crds(2);
02690 NodalCoord.val(1,12)=nd12Crds(0); NodalCoord.val(2,12)=nd12Crds(1); NodalCoord.val(3,12)=nd12Crds(2);
02691 NodalCoord.val(1,13)=nd13Crds(0); NodalCoord.val(2,13)=nd13Crds(1); NodalCoord.val(3,13)=nd13Crds(2);
02692 NodalCoord.val(1,14)=nd14Crds(0); NodalCoord.val(2,14)=nd14Crds(1); NodalCoord.val(3,14)=nd14Crds(2);
02693 NodalCoord.val(1,15)=nd15Crds(0); NodalCoord.val(2,15)=nd15Crds(1); NodalCoord.val(3,15)=nd15Crds(2);
02694 NodalCoord.val(1,16)=nd16Crds(0); NodalCoord.val(2,16)=nd16Crds(1); NodalCoord.val(3,16)=nd16Crds(2);
02695 NodalCoord.val(1,17)=nd17Crds(0); NodalCoord.val(2,17)=nd17Crds(1); NodalCoord.val(3,17)=nd17Crds(2);
02696 NodalCoord.val(1,18)=nd18Crds(0); NodalCoord.val(2,18)=nd18Crds(1); NodalCoord.val(3,18)=nd18Crds(2);
02697 NodalCoord.val(1,19)=nd19Crds(0); NodalCoord.val(2,19)=nd19Crds(1); NodalCoord.val(3,19)=nd19Crds(2);
02698 NodalCoord.val(1,20)=nd20Crds(0); NodalCoord.val(2,20)=nd20Crds(1); NodalCoord.val(3,20)=nd20Crds(2);
02699 NodalCoord.val(1,21)=nd21Crds(0); NodalCoord.val(2,21)=nd21Crds(1); NodalCoord.val(3,21)=nd21Crds(2);
02700 NodalCoord.val(1,22)=nd22Crds(0); NodalCoord.val(2,22)=nd22Crds(1); NodalCoord.val(3,22)=nd22Crds(2);
02701 NodalCoord.val(1,23)=nd23Crds(0); NodalCoord.val(2,23)=nd23Crds(1); NodalCoord.val(3,23)=nd23Crds(2);
02702 NodalCoord.val(1,24)=nd24Crds(0); NodalCoord.val(2,24)=nd24Crds(1); NodalCoord.val(3,24)=nd24Crds(2);
02703 NodalCoord.val(1,25)=nd25Crds(0); NodalCoord.val(2,25)=nd25Crds(1); NodalCoord.val(3,25)=nd25Crds(2);
02704 NodalCoord.val(1,26)=nd26Crds(0); NodalCoord.val(2,26)=nd26Crds(1); NodalCoord.val(3,26)=nd26Crds(2);
02705 NodalCoord.val(1,27)=nd27Crds(0); NodalCoord.val(2,27)=nd27Crds(1); NodalCoord.val(3,27)=nd27Crds(2);
02706
02707
02708
02709 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02710 {
02711 r = get_Gauss_p_c( r_integration_order, GP_c_r );
02712 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02713 {
02714 s = get_Gauss_p_c( s_integration_order, GP_c_s );
02715 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02716 {
02717 t = get_Gauss_p_c( t_integration_order, GP_c_t );
02718
02719
02720 where =
02721 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
02722
02723
02724 H = H_3D(r,s,t);
02725
02726
02727 for (int encount=1 ; encount <= nodes_in_brick; encount++ )
02728
02729 {
02730
02731
02732
02733 matpointCoord.val(1,where+1) +=NodalCoord.val(1,encount) * H.val(encount*3-2,1);
02734
02735 matpointCoord.val(2,where+1) +=NodalCoord.val(2,encount) * H.val(encount*3-1,2);
02736 matpointCoord.val(3,where+1) +=NodalCoord.val(3,encount) * H.val(encount*3-0,3);
02737 }
02738
02739
02740
02741
02742
02743
02744 GaussCoord(where*3+1) = matpointCoord.val(1,where+1);
02745 GaussCoord(where*3+2) = matpointCoord.val(2,where+1);
02746 GaussCoord(where*3+3) = matpointCoord.val(3,where+1);
02747
02748
02749 }
02750 }
02751 }
02752
02753
02754 }
02755
02756
02758
02759
02760
02761
02762
02763 void TwentySevenNodeBrick::reportTensorF(FILE * fp)
02764 {
02765
02766
02767
02768
02769 double r = 0.0;
02770 double s = 0.0;
02771 double t = 0.0;
02772
02773 short where = 0;
02774
02775
02776 static const int dim[] = {3, 27};
02777 tensor NodalCoord(2, dim, 0.0);
02778 tensor matpointCoord(2, dim, 0.0);
02779 int h_dim[] = {81,3};
02780
02781 tensor H(2, h_dim, 0.0);
02782
02783
02784
02785
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804 const Vector &nd1Crds = theNodes[0]->getCrds();
02805 const Vector &nd2Crds = theNodes[1]->getCrds();
02806 const Vector &nd3Crds = theNodes[2]->getCrds();
02807 const Vector &nd4Crds = theNodes[3]->getCrds();
02808 const Vector &nd5Crds = theNodes[4]->getCrds();
02809 const Vector &nd6Crds = theNodes[5]->getCrds();
02810 const Vector &nd7Crds = theNodes[6]->getCrds();
02811 const Vector &nd8Crds = theNodes[7]->getCrds();
02812 const Vector &nd9Crds = theNodes[8]->getCrds();
02813 const Vector &nd10Crds = theNodes[9]->getCrds();
02814 const Vector &nd11Crds = theNodes[10]->getCrds();
02815 const Vector &nd12Crds = theNodes[11]->getCrds();
02816 const Vector &nd13Crds = theNodes[12]->getCrds();
02817 const Vector &nd14Crds = theNodes[13]->getCrds();
02818 const Vector &nd15Crds = theNodes[14]->getCrds();
02819 const Vector &nd16Crds = theNodes[15]->getCrds();
02820 const Vector &nd17Crds = theNodes[16]->getCrds();
02821 const Vector &nd18Crds = theNodes[17]->getCrds();
02822 const Vector &nd19Crds = theNodes[18]->getCrds();
02823 const Vector &nd20Crds = theNodes[19]->getCrds();
02824 const Vector &nd21Crds = theNodes[20]->getCrds();
02825 const Vector &nd22Crds = theNodes[21]->getCrds();
02826 const Vector &nd23Crds = theNodes[22]->getCrds();
02827 const Vector &nd24Crds = theNodes[23]->getCrds();
02828 const Vector &nd25Crds = theNodes[24]->getCrds();
02829 const Vector &nd26Crds = theNodes[25]->getCrds();
02830 const Vector &nd27Crds = theNodes[26]->getCrds();
02831
02832
02833
02834 NodalCoord.val(1,1)=nd1Crds(0); NodalCoord.val(2,1)=nd1Crds(1); NodalCoord.val(3,1)=nd1Crds(2);
02835 NodalCoord.val(1,2)=nd2Crds(0); NodalCoord.val(2,2)=nd2Crds(1); NodalCoord.val(3,2)=nd2Crds(2);
02836 NodalCoord.val(1,3)=nd3Crds(0); NodalCoord.val(2,3)=nd3Crds(1); NodalCoord.val(3,3)=nd3Crds(2);
02837 NodalCoord.val(1,4)=nd4Crds(0); NodalCoord.val(2,4)=nd4Crds(1); NodalCoord.val(3,4)=nd4Crds(2);
02838 NodalCoord.val(1,5)=nd5Crds(0); NodalCoord.val(2,5)=nd5Crds(1); NodalCoord.val(3,5)=nd5Crds(2);
02839 NodalCoord.val(1,6)=nd6Crds(0); NodalCoord.val(2,6)=nd6Crds(1); NodalCoord.val(3,6)=nd6Crds(2);
02840 NodalCoord.val(1,7)=nd7Crds(0); NodalCoord.val(2,7)=nd7Crds(1); NodalCoord.val(3,7)=nd7Crds(2);
02841 NodalCoord.val(1,8)=nd8Crds(0); NodalCoord.val(2,8)=nd8Crds(1); NodalCoord.val(3,8)=nd8Crds(2);
02842 NodalCoord.val(1,9)=nd9Crds(0); NodalCoord.val(2,9)=nd9Crds(1); NodalCoord.val(3,9)=nd9Crds(2);
02843 NodalCoord.val(1,10)=nd10Crds(0); NodalCoord.val(2,10)=nd10Crds(1); NodalCoord.val(3,10)=nd10Crds(2);
02844 NodalCoord.val(1,11)=nd11Crds(0); NodalCoord.val(2,11)=nd11Crds(1); NodalCoord.val(3,11)=nd11Crds(2);
02845 NodalCoord.val(1,12)=nd12Crds(0); NodalCoord.val(2,12)=nd12Crds(1); NodalCoord.val(3,12)=nd12Crds(2);
02846 NodalCoord.val(1,13)=nd13Crds(0); NodalCoord.val(2,13)=nd13Crds(1); NodalCoord.val(3,13)=nd13Crds(2);
02847 NodalCoord.val(1,14)=nd14Crds(0); NodalCoord.val(2,14)=nd14Crds(1); NodalCoord.val(3,14)=nd14Crds(2);
02848 NodalCoord.val(1,15)=nd15Crds(0); NodalCoord.val(2,15)=nd15Crds(1); NodalCoord.val(3,15)=nd15Crds(2);
02849 NodalCoord.val(1,16)=nd16Crds(0); NodalCoord.val(2,16)=nd16Crds(1); NodalCoord.val(3,16)=nd16Crds(2);
02850 NodalCoord.val(1,17)=nd17Crds(0); NodalCoord.val(2,17)=nd17Crds(1); NodalCoord.val(3,17)=nd17Crds(2);
02851 NodalCoord.val(1,18)=nd18Crds(0); NodalCoord.val(2,18)=nd18Crds(1); NodalCoord.val(3,18)=nd18Crds(2);
02852 NodalCoord.val(1,19)=nd19Crds(0); NodalCoord.val(2,19)=nd19Crds(1); NodalCoord.val(3,19)=nd19Crds(2);
02853 NodalCoord.val(1,20)=nd20Crds(0); NodalCoord.val(2,20)=nd20Crds(1); NodalCoord.val(3,20)=nd20Crds(2);
02854 NodalCoord.val(1,21)=nd21Crds(0); NodalCoord.val(2,21)=nd21Crds(1); NodalCoord.val(3,21)=nd21Crds(2);
02855 NodalCoord.val(1,22)=nd22Crds(0); NodalCoord.val(2,22)=nd22Crds(1); NodalCoord.val(3,22)=nd22Crds(2);
02856 NodalCoord.val(1,23)=nd23Crds(0); NodalCoord.val(2,23)=nd23Crds(1); NodalCoord.val(3,23)=nd23Crds(2);
02857 NodalCoord.val(1,24)=nd24Crds(0); NodalCoord.val(2,24)=nd24Crds(1); NodalCoord.val(3,24)=nd24Crds(2);
02858 NodalCoord.val(1,25)=nd25Crds(0); NodalCoord.val(2,25)=nd25Crds(1); NodalCoord.val(3,25)=nd25Crds(2);
02859 NodalCoord.val(1,26)=nd26Crds(0); NodalCoord.val(2,26)=nd26Crds(1); NodalCoord.val(3,26)=nd26Crds(2);
02860 NodalCoord.val(1,27)=nd27Crds(0); NodalCoord.val(2,27)=nd27Crds(1); NodalCoord.val(3,27)=nd27Crds(2);
02861
02862 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
02863 {
02864 r = get_Gauss_p_c( r_integration_order, GP_c_r );
02865 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
02866 {
02867 s = get_Gauss_p_c( s_integration_order, GP_c_s );
02868 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
02869 {
02870 t = get_Gauss_p_c( t_integration_order, GP_c_t );
02871
02872
02873 where =
02874 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
02875
02876
02877 H = H_3D(r,s,t);
02878
02879 for (int encount=1 ; encount <= nodes_in_brick ; encount++ )
02880
02881 {
02882
02883
02884
02885 matpointCoord.val(1,where+1) +=NodalCoord.val(1,encount) * H.val(encount*3-2,1);
02886
02887 matpointCoord.val(2,where+1) +=NodalCoord.val(2,encount) * H.val(encount*3-1,2);
02888 matpointCoord.val(3,where+1) +=NodalCoord.val(3,encount) * H.val(encount*3-0,3);
02889
02890 }
02891
02892 fprintf(fp, "gauss point# %d %+.6e %+.6e %+.6e \n", where+1,
02893 matpointCoord.val(1,where+1),
02894 matpointCoord.val(2,where+1),
02895 matpointCoord.val(3,where+1));
02896
02897
02898
02899
02900 }
02901 }
02902 }
02903
02904 }
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914 int TwentySevenNodeBrick::getNumExternalNodes () const
02915 {
02916 return nodes_in_brick;
02917 }
02918
02919
02920
02921 const ID& TwentySevenNodeBrick::getExternalNodes ()
02922 {
02923 return connectedExternalNodes;
02924 }
02925
02926 Node **
02927 TwentySevenNodeBrick::getNodePtrs(void)
02928 {
02929 return theNodes;
02930 }
02931
02932
02933 int TwentySevenNodeBrick::getNumDOF ()
02934 {
02935 return 3*nodes_in_brick;
02936 }
02937
02938
02939 void TwentySevenNodeBrick::setDomain (Domain *theDomain)
02940 {
02941
02942 if (theDomain == 0) {
02943 theNodes[0] = 0;
02944 theNodes[1] = 0;
02945 theNodes[2] = 0;
02946 theNodes[3] = 0;
02947 theNodes[4] = 0;
02948 theNodes[5] = 0;
02949 theNodes[6] = 0;
02950 theNodes[7] = 0;
02951 theNodes[8] = 0;
02952 theNodes[9] = 0;
02953 theNodes[10] = 0;
02954 theNodes[11] = 0;
02955 theNodes[12] = 0;
02956 theNodes[13] = 0;
02957 theNodes[14] = 0;
02958 theNodes[15] = 0;
02959 theNodes[16] = 0;
02960 theNodes[17] = 0;
02961 theNodes[18] = 0;
02962 theNodes[19] = 0;
02963 theNodes[20] = 0;
02964 theNodes[21] = 0;
02965 theNodes[22] = 0;
02966 theNodes[23] = 0;
02967 theNodes[24] = 0;
02968 theNodes[25] = 0;
02969 theNodes[26] = 0;
02970
02971 }
02972
02973
02974 else {
02975 int Nd1 = connectedExternalNodes(0);
02976 int Nd2 = connectedExternalNodes(1);
02977 int Nd3 = connectedExternalNodes(2);
02978 int Nd4 = connectedExternalNodes(3);
02979
02980
02981 int Nd5 = connectedExternalNodes(4);
02982 int Nd6 = connectedExternalNodes(5);
02983 int Nd7 = connectedExternalNodes(6);
02984 int Nd8 = connectedExternalNodes(7);
02985 int Nd9 = connectedExternalNodes(8);
02986 int Nd10 = connectedExternalNodes(9);
02987 int Nd11 = connectedExternalNodes(10);
02988 int Nd12 = connectedExternalNodes(11);
02989 int Nd13 = connectedExternalNodes(12);
02990 int Nd14 = connectedExternalNodes(13);
02991 int Nd15 = connectedExternalNodes(14);
02992 int Nd16 = connectedExternalNodes(15);
02993 int Nd17 = connectedExternalNodes(16);
02994 int Nd18 = connectedExternalNodes(17);
02995 int Nd19 = connectedExternalNodes(18);
02996 int Nd20 = connectedExternalNodes(19);
02997 int Nd21 = connectedExternalNodes(20);
02998 int Nd22 = connectedExternalNodes(21);
02999 int Nd23 = connectedExternalNodes(22);
03000 int Nd24 = connectedExternalNodes(23);
03001 int Nd25 = connectedExternalNodes(24);
03002 int Nd26 = connectedExternalNodes(25);
03003 int Nd27 = connectedExternalNodes(26);
03004
03005
03006 theNodes[0] = theDomain->getNode(Nd1);
03007 theNodes[1] = theDomain->getNode(Nd2);
03008 theNodes[2] = theDomain->getNode(Nd3);
03009 theNodes[3] = theDomain->getNode(Nd4);
03010 theNodes[4] = theDomain->getNode(Nd5);
03011 theNodes[5] = theDomain->getNode(Nd6);
03012 theNodes[6] = theDomain->getNode(Nd7);
03013 theNodes[7] = theDomain->getNode(Nd8);
03014 theNodes[8] = theDomain->getNode(Nd9);
03015 theNodes[9] = theDomain->getNode(Nd10);
03016 theNodes[10] = theDomain->getNode(Nd11);
03017 theNodes[11] = theDomain->getNode(Nd12);
03018 theNodes[12] = theDomain->getNode(Nd13);
03019 theNodes[13] = theDomain->getNode(Nd14);
03020 theNodes[14] = theDomain->getNode(Nd15);
03021 theNodes[15] = theDomain->getNode(Nd16);
03022 theNodes[16] = theDomain->getNode(Nd17);
03023 theNodes[17] = theDomain->getNode(Nd18);
03024 theNodes[18] = theDomain->getNode(Nd19);
03025 theNodes[19] = theDomain->getNode(Nd20);
03026 theNodes[20] = theDomain->getNode(Nd21);
03027 theNodes[21] = theDomain->getNode(Nd22);
03028 theNodes[22] = theDomain->getNode(Nd23);
03029 theNodes[23] = theDomain->getNode(Nd24);
03030 theNodes[24] = theDomain->getNode(Nd25);
03031 theNodes[25] = theDomain->getNode(Nd26);
03032 theNodes[26] = theDomain->getNode(Nd27);
03033
03034
03035 if (theNodes[0] == 0 || theNodes[1] == 0 || theNodes[2] == 0 || theNodes[3] == 0 ||
03036 theNodes[4] == 0 || theNodes[5] == 0 || theNodes[6] == 0 || theNodes[7] == 0 ||
03037 theNodes[8] == 0 || theNodes[9] == 0 || theNodes[10] == 0 || theNodes[11] == 0 ||
03038 theNodes[12] == 0 || theNodes[13] == 0 || theNodes[14] == 0 || theNodes[15] == 0 ||
03039 theNodes[16] == 0 || theNodes[17] == 0 || theNodes[18] == 0 || theNodes[19] == 0 ||
03040 theNodes[20] == 0 || theNodes[21] == 0 || theNodes[22] == 0 || theNodes[23] == 0 ||
03041 theNodes[24] == 0 || theNodes[25] == 0 || theNodes[26] == 0) {
03042
03043 opserr << "FATAL ERROR TwentySevenNodeBrick (tag: " << this->getTag() << " ), node not found in domain\n";
03044 exit(-1);
03045 }
03046
03047 int dofNd1 = theNodes[0]->getNumberDOF();
03048 int dofNd2 = theNodes[1]->getNumberDOF();
03049 int dofNd3 = theNodes[2]->getNumberDOF();
03050 int dofNd4 = theNodes[3]->getNumberDOF();
03051 int dofNd5 = theNodes[4]->getNumberDOF();
03052 int dofNd6 = theNodes[5]->getNumberDOF();
03053 int dofNd7 = theNodes[6]->getNumberDOF();
03054 int dofNd8 = theNodes[7]->getNumberDOF();
03055 int dofNd9 = theNodes[8]->getNumberDOF();
03056 int dofNd10 = theNodes[9]->getNumberDOF();
03057 int dofNd11 = theNodes[10]->getNumberDOF();
03058 int dofNd12 = theNodes[11]->getNumberDOF();
03059 int dofNd13 = theNodes[12]->getNumberDOF();
03060 int dofNd14 = theNodes[13]->getNumberDOF();
03061 int dofNd15 = theNodes[14]->getNumberDOF();
03062 int dofNd16 = theNodes[15]->getNumberDOF();
03063 int dofNd17 = theNodes[16]->getNumberDOF();
03064 int dofNd18 = theNodes[17]->getNumberDOF();
03065 int dofNd19 = theNodes[18]->getNumberDOF();
03066 int dofNd20 = theNodes[19]->getNumberDOF();
03067 int dofNd21 = theNodes[20]->getNumberDOF();
03068 int dofNd22 = theNodes[21]->getNumberDOF();
03069 int dofNd23 = theNodes[22]->getNumberDOF();
03070 int dofNd24 = theNodes[23]->getNumberDOF();
03071 int dofNd25 = theNodes[24]->getNumberDOF();
03072 int dofNd26 = theNodes[25]->getNumberDOF();
03073 int dofNd27 = theNodes[26]->getNumberDOF();
03074
03075
03076 if (dofNd1 != 3 || dofNd2 != 3 || dofNd3 != 3 || dofNd4 != 3 ||
03077 dofNd5 != 3 || dofNd6 != 3 || dofNd7 != 3 || dofNd8 != 3 ||
03078 dofNd9 != 3 || dofNd10 != 3 || dofNd11 != 3 || dofNd12 != 3 ||
03079 dofNd13 != 3 || dofNd14 != 3 || dofNd15 != 3 || dofNd16 != 3 ||
03080 dofNd17 != 3 || dofNd18 != 3 || dofNd19 != 3 || dofNd20 != 3 ||
03081 dofNd21 != 3 || dofNd22 != 3 || dofNd23 != 3 || dofNd24 != 3 ||
03082 dofNd25 != 3 || dofNd26 != 3 || dofNd27 != 3){
03083 opserr << "FATAL ERROR TwentySevenNodeBrick (tag: " << this->getTag() <<
03084 "), has differing number of DOFs at its nodes\n";
03085 exit(-1);
03086 }
03087 this->DomainComponent::setDomain(theDomain);
03088 }
03089 }
03090
03091
03092 int TwentySevenNodeBrick::commitState ()
03093 {
03094 int retVal = 0;
03095
03096
03097 if ((retVal = this->Element::commitState()) != 0) {
03098 opserr << "TwentySevenNodeBrick::commitState () - failed in base class";
03099 }
03100
03101
03102
03103 int i;
03104
03105
03106
03107 int count = r_integration_order* s_integration_order * t_integration_order;
03108
03109
03110
03111
03112
03113
03114 Vector pp = getResistingForce();
03115
03116
03117
03118 for (i = 0; i < count; i++)
03119
03120 {
03121 retVal += matpoint[i]->commitState();
03122
03123 stresstensor st;
03124 stresstensor prin;
03125 straintensor stn;
03126 straintensor stnprin;
03127
03128 st = matpoint[i]->getStressTensor();
03129 prin = st.principal();
03130 stn = matpoint[i]->getStrainTensor();
03131 stnprin = stn.principal();
03132
03133
03134
03135
03136
03137
03138
03139
03140
03141
03142 double p = -1*( prin.cval(1, 1)+ prin.cval(2, 2) +prin.cval(3, 3) )/3.0;
03143 double ev = -1*( stnprin.cval(1, 1)+ stnprin.cval(2, 2) + stnprin.cval(3, 3) )/3.0;
03144
03145
03146
03147
03148
03149
03150 double q;
03151
03152 if ( fabs(prin.cval(1, 1) - prin.cval(2, 2) ) <= 0.001 )
03153 {
03154 q = prin.cval(1, 1) - prin.cval(3, 3);
03155
03156 }
03157 else
03158 q = prin.cval(3, 3) - prin.cval(1, 1);
03159
03160
03161
03162
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181
03182
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205
03206
03207
03208
03209 }
03210
03211
03212
03213
03214
03215
03216
03217
03218 return retVal;
03219 }
03220
03221
03222 int TwentySevenNodeBrick::revertToLastCommit ()
03223 {
03224
03225 int i;
03226
03227 int retVal = 0;
03228
03229
03230 int count = r_integration_order* s_integration_order * t_integration_order;
03231
03232
03233
03234
03235
03236
03237 for (i = 0; i < count; i++)
03238 retVal += matpoint[i]->revertToLastCommit();
03239
03240
03241 return retVal;
03242 }
03243
03244
03245 int TwentySevenNodeBrick::revertToStart ()
03246 {
03247 int i;
03248 int retVal = 0;
03249
03250
03251
03252
03253
03254
03255
03256
03257 int count = r_integration_order* s_integration_order * t_integration_order;
03258
03259 for (i = 0; i < count; i++)
03260 retVal += matpoint[i]->revertToStart();
03261
03262
03263 return retVal;
03264
03265
03266 }
03267
03268
03269
03270 const Matrix &TwentySevenNodeBrick::getTangentStiff ()
03271 {
03272 tensor stifftensor = getStiffnessTensor();
03273 int Ki=0;
03274 int Kj=0;
03275
03276 for ( int i=1 ; i<=nodes_in_brick ; i++ )
03277 {
03278 for ( int j=1 ; j<=nodes_in_brick ; j++ )
03279 {
03280 for ( int k=1 ; k<=3 ; k++ )
03281 {
03282 for ( int l=1 ; l<=3 ; l++ )
03283 {
03284 Ki = k+3*(i-1);
03285 Kj = l+3*(j-1);
03286 K( Ki-1 , Kj-1 ) = stifftensor.cval(i,k,l,j);
03287 }
03288 }
03289 }
03290 }
03291
03292
03293
03294 return K;
03295 }
03296
03297
03298 const Matrix &TwentySevenNodeBrick::getInitialStiff ()
03299 {
03300
03301 return this->getTangentStiff();
03302 }
03303
03304
03305
03306
03307 const Matrix &TwentySevenNodeBrick::getConsMass ()
03308 {
03309 tensor masstensor = getMassTensor();
03310
03311
03312
03313
03314
03315 double column_mass;
03316
03317 for ( int i=1 ; i<=nodes_in_brick*3 ; i++ )
03318 {
03319 column_mass = 0.0;
03320 for ( int j=1 ; j<=nodes_in_brick*3 ; j++ )
03321 {
03322
03323
03324
03325 column_mass += masstensor.cval(i,j);
03326 M( i-1 , j-1 ) = 0;
03327
03328
03329
03330 }
03331 M( i-1 , i-1 ) = column_mass;
03332
03333 }
03334
03335
03336
03337
03338
03339 return M;
03340 }
03341
03342
03343
03344
03345 const Matrix &TwentySevenNodeBrick::getMass ()
03346 {
03347 tensor masstensor = getMassTensor();
03348
03349
03350
03351
03352
03353
03354
03355 for ( int i=1 ; i<=nodes_in_brick*3 ; i++ )
03356 {
03357
03358 for ( int j=1 ; j<=nodes_in_brick*3 ; j++ )
03359 {
03360 M( i-1 , j-1 ) = masstensor.cval(i,j);
03361
03362
03363
03364
03365
03366
03367 }
03368
03369
03370 }
03371
03372
03373
03374
03375
03376 return M;
03377 }
03378
03379
03380 void TwentySevenNodeBrick::zeroLoad(void)
03381 {
03382 Q.Zero();
03383 }
03384
03385
03386
03387 int
03388 TwentySevenNodeBrick::addLoad(ElementalLoad *theLoad, double loadFactor)
03389 {
03390
03391
03392
03393 int type;
03394 const Vector &data = theLoad->getData(type, loadFactor);
03395 if (type == LOAD_TAG_BrickSelfWeight) {
03396
03397 Vector bforce(81);
03398
03399
03400 if (rho == 0.0)
03401 return 0;
03402
03403 Vector ba(81), bfx(3);
03404 bfx(0) = bf(0) * loadFactor;
03405 bfx(1) = bf(1) * loadFactor;
03406 bfx(2) = bf(2) * loadFactor;
03407
03408 ba( 0) = bfx(0);
03409 ba( 1) = bfx(1);
03410 ba( 2) = bfx(2);
03411 ba( 3) = bfx(0);
03412 ba( 4) = bfx(1);
03413 ba( 5) = bfx(2);
03414 ba( 6) = bfx(0);
03415 ba( 7) = bfx(1);
03416 ba( 8) = bfx(2);
03417 ba( 9) = bfx(0);
03418 ba(10) = bfx(1);
03419 ba(11) = bfx(2);
03420 ba(12) = bfx(0);
03421 ba(13) = bfx(1);
03422 ba(14) = bfx(2);
03423 ba(15) = bfx(0);
03424 ba(16) = bfx(1);
03425 ba(17) = bfx(2);
03426 ba(18) = bfx(0);
03427 ba(19) = bfx(1);
03428 ba(20) = bfx(2);
03429 ba(21) = bfx(0);
03430 ba(22) = bfx(1);
03431 ba(23) = bfx(2);
03432 ba(24) = bfx(0);
03433 ba(25) = bfx(1);
03434 ba(26) = bfx(2);
03435 ba(27) = bfx(0);
03436 ba(28) = bfx(1);
03437 ba(29) = bfx(2);
03438 ba(30) = bfx(0);
03439 ba(31) = bfx(1);
03440 ba(32) = bfx(2);
03441 ba(33) = bfx(0);
03442 ba(34) = bfx(1);
03443 ba(35) = bfx(2);
03444 ba(36) = bfx(0);
03445 ba(37) = bfx(1);
03446 ba(38) = bfx(2);
03447 ba(39) = bfx(0);
03448 ba(40) = bfx(1);
03449 ba(41) = bfx(2);
03450 ba(42) = bfx(0);
03451 ba(43) = bfx(1);
03452 ba(44) = bfx(2);
03453 ba(45) = bfx(0);
03454 ba(46) = bfx(1);
03455 ba(47) = bfx(2);
03456 ba(48) = bfx(0);
03457 ba(49) = bfx(1);
03458 ba(50) = bfx(2);
03459 ba(51) = bfx(0);
03460 ba(52) = bfx(1);
03461 ba(53) = bfx(2);
03462 ba(54) = bfx(0);
03463 ba(55) = bfx(1);
03464 ba(56) = bfx(2);
03465 ba(57) = bfx(0);
03466 ba(58) = bfx(1);
03467 ba(59) = bfx(2);
03468 ba(60) = bfx(0);
03469 ba(61) = bfx(1);
03470 ba(62) = bfx(2);
03471 ba(63) = bfx(0);
03472 ba(64) = bfx(1);
03473 ba(65) = bfx(2);
03474 ba(66) = bfx(0);
03475 ba(67) = bfx(1);
03476 ba(68) = bfx(2);
03477 ba(69) = bfx(0);
03478 ba(70) = bfx(1);
03479 ba(71) = bfx(2);
03480 ba(72) = bfx(0);
03481 ba(73) = bfx(1);
03482 ba(74) = bfx(2);
03483 ba(75) = bfx(0);
03484 ba(76) = bfx(1);
03485 ba(77) = bfx(2);
03486 ba(78) = bfx(0);
03487 ba(79) = bfx(1);
03488 ba(80) = bfx(2);
03489
03490
03491
03492
03493
03494 this->getMass();
03495 bforce.addMatrixVector(0.0, M, ba, 1.0);
03496 Q.addVector(1.0, bforce, 1.0);
03497 } else {
03498 opserr << "TwentySevenNodeBrick::addLoad() - 20NodeBrick " << this->getTag() <<
03499 ",load type " << type << " unknown\n";
03500 return -1;
03501 }
03502
03503 return 0;
03504 }
03505
03506
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522
03523 int TwentySevenNodeBrick::addInertiaLoadToUnbalance(const Vector &accel)
03524 {
03525
03526 if (rho == 0.0)
03527 return 0;
03528
03529
03530 const Vector &Raccel1 = theNodes[0]->getRV(accel);
03531 const Vector &Raccel2 = theNodes[1]->getRV(accel);
03532 const Vector &Raccel3 = theNodes[2]->getRV(accel);
03533 const Vector &Raccel4 = theNodes[3]->getRV(accel);
03534 const Vector &Raccel5 = theNodes[4]->getRV(accel);
03535 const Vector &Raccel6 = theNodes[5]->getRV(accel);
03536 const Vector &Raccel7 = theNodes[6]->getRV(accel);
03537 const Vector &Raccel8 = theNodes[7]->getRV(accel);
03538 const Vector &Raccel9 = theNodes[8]->getRV(accel);
03539 const Vector &Raccel10 = theNodes[9]->getRV(accel);
03540 const Vector &Raccel11 = theNodes[10]->getRV(accel);
03541 const Vector &Raccel12 = theNodes[11]->getRV(accel);
03542 const Vector &Raccel13 = theNodes[12]->getRV(accel);
03543 const Vector &Raccel14 = theNodes[13]->getRV(accel);
03544 const Vector &Raccel15 = theNodes[14]->getRV(accel);
03545 const Vector &Raccel16 = theNodes[15]->getRV(accel);
03546 const Vector &Raccel17 = theNodes[16]->getRV(accel);
03547 const Vector &Raccel18 = theNodes[17]->getRV(accel);
03548 const Vector &Raccel19 = theNodes[18]->getRV(accel);
03549 const Vector &Raccel20 = theNodes[19]->getRV(accel);
03550 const Vector &Raccel21 = theNodes[20]->getRV(accel);
03551 const Vector &Raccel22 = theNodes[21]->getRV(accel);
03552 const Vector &Raccel23 = theNodes[22]->getRV(accel);
03553 const Vector &Raccel24 = theNodes[23]->getRV(accel);
03554 const Vector &Raccel25 = theNodes[24]->getRV(accel);
03555 const Vector &Raccel26 = theNodes[25]->getRV(accel);
03556 const Vector &Raccel27 = theNodes[26]->getRV(accel);
03557
03558
03559 if (3 != Raccel1.Size() || 3 != Raccel2.Size() || 3 != Raccel3.Size() || 3 != Raccel4.Size() ||
03560 3 != Raccel5.Size() || 3 != Raccel6.Size() || 3 != Raccel7.Size() || 3 != Raccel8.Size() ||
03561 3 != Raccel9.Size() || 3 != Raccel10.Size() || 3 != Raccel11.Size() || 3 != Raccel12.Size()||
03562 3 != Raccel13.Size() || 3 != Raccel14.Size() || 3 != Raccel15.Size() || 3 != Raccel16.Size()||
03563 3 != Raccel17.Size() || 3 != Raccel18.Size() || 3 != Raccel19.Size() || 3 != Raccel20.Size()||
03564 3 != Raccel21.Size() || 3 != Raccel22.Size() || 3 != Raccel23.Size() || 3 != Raccel24.Size()||
03565 3 != Raccel25.Size() || 3 != Raccel26.Size() || 3 != Raccel27.Size()){
03566
03567 opserr << "TwentySevenNodeBrick::addInertiaLoadToUnbalance " <<
03568 "matrix and vector sizes are incompatable\n";
03569 return -1;
03570 }
03571
03572 static Vector ra(81);
03573
03574 ra( 0) = Raccel1(0);
03575 ra( 1) = Raccel1(1);
03576 ra( 2) = Raccel1(2);
03577 ra( 3) = Raccel2(0);
03578 ra( 4) = Raccel2(1);
03579 ra( 5) = Raccel2(2);
03580 ra( 6) = Raccel3(0);
03581 ra( 7) = Raccel3(1);
03582 ra( 8) = Raccel3(2);
03583 ra( 9) = Raccel4(0);
03584 ra(10) = Raccel4(1);
03585 ra(11) = Raccel4(2);
03586 ra(12) = Raccel5(0);
03587 ra(13) = Raccel5(1);
03588 ra(14) = Raccel5(2);
03589 ra(15) = Raccel6(0);
03590 ra(16) = Raccel6(1);
03591 ra(17) = Raccel6(2);
03592 ra(18) = Raccel7(0);
03593 ra(19) = Raccel7(1);
03594 ra(20) = Raccel7(2);
03595 ra(21) = Raccel8(0);
03596 ra(22) = Raccel8(1);
03597 ra(23) = Raccel8(2);
03598 ra(24) = Raccel9(0);
03599 ra(25) = Raccel9(1);
03600 ra(26) = Raccel9(2);
03601 ra(27) = Raccel10(0);
03602 ra(28) = Raccel10(1);
03603 ra(29) = Raccel10(2);
03604 ra(30) = Raccel11(0);
03605 ra(31) = Raccel11(1);
03606 ra(32) = Raccel11(2);
03607 ra(33) = Raccel12(0);
03608 ra(34) = Raccel12(1);
03609 ra(35) = Raccel12(2);
03610 ra(36) = Raccel13(0);
03611 ra(37) = Raccel13(1);
03612 ra(38) = Raccel13(2);
03613 ra(39) = Raccel14(0);
03614 ra(40) = Raccel14(1);
03615 ra(41) = Raccel14(2);
03616 ra(42) = Raccel15(0);
03617 ra(43) = Raccel15(1);
03618 ra(44) = Raccel15(2);
03619 ra(45) = Raccel16(0);
03620 ra(46) = Raccel16(1);
03621 ra(47) = Raccel16(2);
03622 ra(48) = Raccel17(0);
03623 ra(49) = Raccel17(1);
03624 ra(50) = Raccel17(2);
03625 ra(51) = Raccel18(0);
03626 ra(52) = Raccel18(1);
03627 ra(53) = Raccel18(2);
03628 ra(54) = Raccel19(0);
03629 ra(55) = Raccel19(1);
03630 ra(56) = Raccel19(2);
03631 ra(57) = Raccel20(0);
03632 ra(58) = Raccel20(1);
03633 ra(59) = Raccel20(2);
03634 ra(60) = Raccel21(0);
03635 ra(61) = Raccel21(1);
03636 ra(62) = Raccel21(2);
03637 ra(63) = Raccel22(0);
03638 ra(64) = Raccel22(1);
03639 ra(65) = Raccel22(2);
03640 ra(66) = Raccel23(0);
03641 ra(67) = Raccel23(1);
03642 ra(68) = Raccel23(2);
03643 ra(69) = Raccel24(0);
03644 ra(70) = Raccel24(1);
03645 ra(71) = Raccel24(2);
03646 ra(72) = Raccel25(0);
03647 ra(73) = Raccel25(1);
03648 ra(74) = Raccel25(2);
03649 ra(75) = Raccel26(0);
03650 ra(76) = Raccel26(1);
03651 ra(77) = Raccel26(2);
03652 ra(78) = Raccel27(0);
03653 ra(79) = Raccel27(1);
03654 ra(80) = Raccel27(2);
03655
03656
03657
03658
03659
03660
03661
03662
03663
03664
03665
03666
03667
03668
03669
03670
03671 Q.addMatrixVector(1.0, M, ra, -1.0);
03672
03673 return 0;
03674 }
03675
03676
03677 const Vector TwentySevenNodeBrick::FormEquiBodyForce(void)
03678 {
03679 Vector bforce(81);
03680
03681
03682
03683 if (rho == 0.0)
03684 return bforce;
03685
03686 Vector ba(81);
03687
03688 ba( 0) = bf(0);
03689 ba( 1) = bf(1);
03690 ba( 2) = bf(2);
03691 ba( 3) = bf(0);
03692 ba( 4) = bf(1);
03693 ba( 5) = bf(2);
03694 ba( 6) = bf(0);
03695 ba( 7) = bf(1);
03696 ba( 8) = bf(2);
03697 ba( 9) = bf(0);
03698 ba(10) = bf(1);
03699 ba(11) = bf(2);
03700 ba(12) = bf(0);
03701 ba(13) = bf(1);
03702 ba(14) = bf(2);
03703 ba(15) = bf(0);
03704 ba(16) = bf(1);
03705 ba(17) = bf(2);
03706 ba(18) = bf(0);
03707 ba(19) = bf(1);
03708 ba(20) = bf(2);
03709 ba(21) = bf(0);
03710 ba(22) = bf(1);
03711 ba(23) = bf(2);
03712 ba(24) = bf(0);
03713 ba(25) = bf(1);
03714 ba(26) = bf(2);
03715 ba(27) = bf(0);
03716 ba(28) = bf(1);
03717 ba(29) = bf(2);
03718 ba(30) = bf(0);
03719 ba(31) = bf(1);
03720 ba(32) = bf(2);
03721 ba(33) = bf(0);
03722 ba(34) = bf(1);
03723 ba(35) = bf(2);
03724 ba(36) = bf(0);
03725 ba(37) = bf(1);
03726 ba(38) = bf(2);
03727 ba(39) = bf(0);
03728 ba(40) = bf(1);
03729 ba(41) = bf(2);
03730 ba(42) = bf(0);
03731 ba(43) = bf(1);
03732 ba(44) = bf(2);
03733 ba(45) = bf(0);
03734 ba(46) = bf(1);
03735 ba(47) = bf(2);
03736 ba(48) = bf(0);
03737 ba(49) = bf(1);
03738 ba(50) = bf(2);
03739 ba(51) = bf(0);
03740 ba(52) = bf(1);
03741 ba(53) = bf(2);
03742 ba(54) = bf(0);
03743 ba(55) = bf(1);
03744 ba(56) = bf(2);
03745 ba(57) = bf(0);
03746 ba(58) = bf(1);
03747 ba(59) = bf(2);
03748 ba(60) = bf(0);
03749 ba(61) = bf(1);
03750 ba(62) = bf(2);
03751 ba(63) = bf(0);
03752 ba(64) = bf(1);
03753 ba(65) = bf(2);
03754 ba(66) = bf(0);
03755 ba(67) = bf(1);
03756 ba(68) = bf(2);
03757 ba(69) = bf(0);
03758 ba(70) = bf(1);
03759 ba(71) = bf(2);
03760 ba(72) = bf(0);
03761 ba(73) = bf(1);
03762 ba(74) = bf(2);
03763 ba(75) = bf(0);
03764 ba(76) = bf(1);
03765 ba(77) = bf(2);
03766 ba(78) = bf(0);
03767 ba(79) = bf(1);
03768 ba(80) = bf(2);
03769
03770
03771
03772 bforce.addMatrixVector(0.0, M, ba, 1.0);
03773
03774
03775
03776
03777
03778 return bforce;
03779 }
03780
03781
03782
03783
03784
03785
03786
03787
03788
03789
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802
03803
03804
03805
03806
03807
03808
03809
03810
03811
03812
03813
03814
03815
03816
03817
03818
03819
03820
03821
03822
03823
03824
03825
03826
03827
03828
03829
03830 const Vector &TwentySevenNodeBrick::getResistingForce ()
03831 {
03832 int force_dim[] = {27,3};
03833 tensor nodalforces(2,force_dim,0.0);
03834
03835 nodalforces = nodal_forces();
03836
03837
03838 for (int i = 0; i< nodes_in_brick; i++)
03839 for (int j = 0; j < 3; j++)
03840 P(i *3 + j) = nodalforces.cval(i+1, j+1);
03841
03842
03843
03844
03845 P = P - Q;
03846
03847
03848 return P;
03849 }
03850
03851
03852 const Vector &TwentySevenNodeBrick::getResistingForceIncInertia ()
03853 {
03854
03855 this->getResistingForce();
03856
03857
03858
03859
03860
03861
03862 if (rho != 0.0) {
03863
03864
03865 this->getMass();
03866
03867 const Vector &accel1 = theNodes[0]->getTrialAccel();
03868 const Vector &accel2 = theNodes[1]->getTrialAccel();
03869 const Vector &accel3 = theNodes[2]->getTrialAccel();
03870 const Vector &accel4 = theNodes[3]->getTrialAccel();
03871 const Vector &accel5 = theNodes[4]->getTrialAccel();
03872 const Vector &accel6 = theNodes[5]->getTrialAccel();
03873 const Vector &accel7 = theNodes[6]->getTrialAccel();
03874 const Vector &accel8 = theNodes[7]->getTrialAccel();
03875 const Vector &accel9 = theNodes[8]->getTrialAccel();
03876 const Vector &accel10 = theNodes[9]->getTrialAccel();
03877 const Vector &accel11 = theNodes[10]->getTrialAccel();
03878 const Vector &accel12 = theNodes[11]->getTrialAccel();
03879 const Vector &accel13 = theNodes[12]->getTrialAccel();
03880 const Vector &accel14 = theNodes[13]->getTrialAccel();
03881 const Vector &accel15 = theNodes[14]->getTrialAccel();
03882 const Vector &accel16 = theNodes[15]->getTrialAccel();
03883 const Vector &accel17 = theNodes[16]->getTrialAccel();
03884 const Vector &accel18 = theNodes[17]->getTrialAccel();
03885 const Vector &accel19 = theNodes[18]->getTrialAccel();
03886 const Vector &accel20 = theNodes[19]->getTrialAccel();
03887 const Vector &accel21 = theNodes[20]->getTrialAccel();
03888 const Vector &accel22 = theNodes[21]->getTrialAccel();
03889 const Vector &accel23 = theNodes[22]->getTrialAccel();
03890 const Vector &accel24 = theNodes[23]->getTrialAccel();
03891 const Vector &accel25 = theNodes[24]->getTrialAccel();
03892 const Vector &accel26 = theNodes[25]->getTrialAccel();
03893 const Vector &accel27 = theNodes[26]->getTrialAccel();
03894
03895 static Vector a(81);
03896
03897 a( 0) = accel1(0);
03898 a( 1) = accel1(1);
03899 a( 2) = accel1(2);
03900 a( 3) = accel2(0);
03901 a( 4) = accel2(1);
03902 a( 5) = accel2(2);
03903 a( 6) = accel3(0);
03904 a( 7) = accel3(1);
03905 a( 8) = accel3(2);
03906 a( 9) = accel4(0);
03907 a(10) = accel4(1);
03908 a(11) = accel4(2);
03909 a(12) = accel5(0);
03910 a(13) = accel5(1);
03911 a(14) = accel5(2);
03912 a(15) = accel6(0);
03913 a(16) = accel6(1);
03914 a(17) = accel6(2);
03915 a(18) = accel7(0);
03916 a(19) = accel7(1);
03917 a(20) = accel7(2);
03918 a(21) = accel8(0);
03919 a(22) = accel8(1);
03920 a(23) = accel8(2);
03921 a(24) = accel9(0);
03922 a(25) = accel9(1);
03923 a(26) = accel9(2);
03924 a(27) = accel10(0);
03925 a(28) = accel10(1);
03926 a(29) = accel10(2);
03927 a(30) = accel11(0);
03928 a(31) = accel11(1);
03929 a(32) = accel11(2);
03930 a(33) = accel12(0);
03931 a(34) = accel12(1);
03932 a(35) = accel12(2);
03933 a(36) = accel13(0);
03934 a(37) = accel13(1);
03935 a(38) = accel13(2);
03936 a(39) = accel14(0);
03937 a(40) = accel14(1);
03938 a(41) = accel14(2);
03939 a(42) = accel15(0);
03940 a(43) = accel15(1);
03941 a(44) = accel15(2);
03942 a(45) = accel16(0);
03943 a(46) = accel16(1);
03944 a(47) = accel16(2);
03945 a(48) = accel17(0);
03946 a(49) = accel17(1);
03947 a(50) = accel17(2);
03948 a(51) = accel18(0);
03949 a(52) = accel18(1);
03950 a(53) = accel18(2);
03951 a(54) = accel19(0);
03952 a(55) = accel19(1);
03953 a(56) = accel19(2);
03954 a(57) = accel20(0);
03955 a(58) = accel20(1);
03956 a(59) = accel20(2);
03957 a(60) = accel21(0);
03958 a(61) = accel21(1);
03959 a(62) = accel21(2);
03960 a(63) = accel22(0);
03961 a(64) = accel22(1);
03962 a(65) = accel22(2);
03963 a(66) = accel23(0);
03964 a(67) = accel23(1);
03965 a(68) = accel23(2);
03966 a(69) = accel24(0);
03967 a(70) = accel24(1);
03968 a(71) = accel24(2);
03969 a(72) = accel25(0);
03970 a(73) = accel25(1);
03971 a(74) = accel25(2);
03972 a(75) = accel26(0);
03973 a(76) = accel26(1);
03974 a(77) = accel26(2);
03975 a(78) = accel27(0);
03976 a(79) = accel27(1);
03977 a(80) = accel27(2);
03978
03979
03980
03981
03982 P.addMatrixVector(1.0, M, a, 1.0);
03983
03984
03985 if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
03986 P += this->getRayleighDampingForces();
03987
03988 } else {
03989
03990
03991 if (betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
03992 P += this->getRayleighDampingForces();
03993
03994 }
03995
03996 return P;
03997 }
03998
03999
04000 int TwentySevenNodeBrick::sendSelf (int commitTag, Channel &theChannel)
04001 {
04002
04003 return 0;
04004 }
04005
04006
04007 int TwentySevenNodeBrick::recvSelf (int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
04008 {
04009
04010 return 0;
04011 }
04012
04013
04014
04015 int TwentySevenNodeBrick::displaySelf (Renderer &theViewer, int displayMode, float fact)
04016 {
04017
04018
04019
04020
04021
04022 const Vector &end1Crd = theNodes[0]->getCrds();
04023 const Vector &end2Crd = theNodes[1]->getCrds();
04024 const Vector &end3Crd = theNodes[2]->getCrds();
04025 const Vector &end4Crd = theNodes[3]->getCrds();
04026 const Vector &end5Crd = theNodes[4]->getCrds();
04027 const Vector &end6Crd = theNodes[5]->getCrds();
04028 const Vector &end7Crd = theNodes[6]->getCrds();
04029 const Vector &end8Crd = theNodes[7]->getCrds();
04030
04031 const Vector &end1Disp = theNodes[0]->getDisp();
04032 const Vector &end2Disp = theNodes[1]->getDisp();
04033 const Vector &end3Disp = theNodes[2]->getDisp();
04034 const Vector &end4Disp = theNodes[3]->getDisp();
04035 const Vector &end5Disp = theNodes[4]->getDisp();
04036 const Vector &end6Disp = theNodes[5]->getDisp();
04037 const Vector &end7Disp = theNodes[6]->getDisp();
04038 const Vector &end8Disp = theNodes[7]->getDisp();
04039
04040 static Vector v1(3);
04041 static Vector v2(3);
04042 static Vector v3(3);
04043 static Vector v4(3);
04044 static Vector v5(3);
04045 static Vector v6(3);
04046 static Vector v7(3);
04047 static Vector v8(3);
04048
04049 for (int i = 0; i < 2; i++)
04050 {
04051 v1(i) = end1Crd(i) + end1Disp(i)*fact;
04052 v2(i) = end2Crd(i) + end2Disp(i)*fact;
04053 v3(i) = end3Crd(i) + end3Disp(i)*fact;
04054 v4(i) = end4Crd(i) + end4Disp(i)*fact;
04055 v5(i) = end5Crd(i) + end5Disp(i)*fact;
04056 v6(i) = end6Crd(i) + end6Disp(i)*fact;
04057 v7(i) = end7Crd(i) + end7Disp(i)*fact;
04058 v8(i) = end8Crd(i) + end8Disp(i)*fact;
04059 }
04060
04061 int error = 0;
04062
04063 error += theViewer.drawLine (v1, v2, 1.0, 1.0);
04064 error += theViewer.drawLine (v2, v3, 1.0, 1.0);
04065 error += theViewer.drawLine (v3, v4, 1.0, 1.0);
04066 error += theViewer.drawLine (v4, v1, 1.0, 1.0);
04067
04068 error += theViewer.drawLine (v5, v6, 1.0, 1.0);
04069 error += theViewer.drawLine (v6, v7, 1.0, 1.0);
04070 error += theViewer.drawLine (v7, v8, 1.0, 1.0);
04071 error += theViewer.drawLine (v8, v5, 1.0, 1.0);
04072
04073 error += theViewer.drawLine (v1, v5, 1.0, 1.0);
04074 error += theViewer.drawLine (v2, v6, 1.0, 1.0);
04075 error += theViewer.drawLine (v3, v7, 1.0, 1.0);
04076 error += theViewer.drawLine (v4, v8, 1.0, 1.0);
04077
04078 return error;
04079
04080 }
04081
04082
04083 void TwentySevenNodeBrick::Print(OPS_Stream &s, int flag)
04084 {
04085
04086 s << "TwentySevenNodeBrick, element id: " << this->getTag() << endln;
04087 s << "Connected external nodes: " << connectedExternalNodes;
04088
04089 int total_number_of_Gauss_points = r_integration_order*
04090 s_integration_order*
04091 t_integration_order;
04092 if ( total_number_of_Gauss_points != 0 )
04093 {
04094 theNodes[0]->Print(opserr);
04095 theNodes[1]->Print(opserr);
04096 theNodes[2]->Print(opserr);
04097 theNodes[3]->Print(opserr);
04098 theNodes[4]->Print(opserr);
04099 theNodes[5]->Print(opserr);
04100 theNodes[6]->Print(opserr);
04101 theNodes[7]->Print(opserr);
04102 theNodes[8]->Print(opserr);
04103 theNodes[9]->Print(opserr);
04104 theNodes[10]->Print(opserr);
04105 theNodes[11]->Print(opserr);
04106 theNodes[12]->Print(opserr);
04107 theNodes[13]->Print(opserr);
04108 theNodes[14]->Print(opserr);
04109 theNodes[15]->Print(opserr);
04110 theNodes[16]->Print(opserr);
04111 theNodes[17]->Print(opserr);
04112 theNodes[18]->Print(opserr);
04113 theNodes[19]->Print(opserr);
04114 theNodes[20]->Print(opserr);
04115 theNodes[21]->Print(opserr);
04116 theNodes[22]->Print(opserr);
04117 theNodes[23]->Print(opserr);
04118 theNodes[24]->Print(opserr);
04119 theNodes[25]->Print(opserr);
04120 theNodes[26]->Print(opserr);
04121
04122
04123 }
04124 s << "Element mass density: " << rho << endln << endln;
04125 s << "Material model: " << endln;
04126
04127 for( int GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
04128 {
04129 for( int GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
04130 {
04131 for( int GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
04132 {
04133
04134
04135 short where =
04136 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
04137
04138 s << "\n where = " << where << endln;
04139 s << " GP_c_r= " << GP_c_r << "GP_c_s = " << GP_c_s << " GP_c_t = " << GP_c_t << endln;
04140 matpoint[where]->report("Material Point\n");
04141
04142
04143
04144 }
04145 }
04146 }
04147
04148 }
04149
04150
04151 Response * TwentySevenNodeBrick::setResponse (const char **argv, int argc, Information &eleInformation)
04152 {
04153
04154 if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0)
04155 return new ElementResponse(this, 1, P);
04156
04157
04158 else if (strcmp(argv[0],"stiff") == 0 || strcmp(argv[0],"stiffness") == 0)
04159 return new ElementResponse(this, 5, K);
04160
04161
04162 else if (strcmp(argv[0],"plastic") == 0 || strcmp(argv[0],"plastified") == 0)
04163 {
04165
04166
04167
04168
04169
04170
04171
04172
04173
04174
04175
04176
04177
04178
04179 return new ElementResponse(this, 2, InfoPlastic);
04180 }
04181
04182
04183 else if (strcmp(argv[0],"PileM") == 0 || strcmp(argv[0],"PileM") == 0)
04184 {
04185 return new ElementResponse(this, 3, InfoStress);
04186 }
04187
04188 else if (strcmp(argv[0],"stress") == 0 || strcmp(argv[0],"stresses") == 0)
04189 {
04190 return new ElementResponse(this, 4, InfoStress);
04191 }
04192
04193
04194 else if (strcmp(argv[0],"pq") == 0 || strcmp(argv[0],"PQ") == 0)
04195 {
04196 return new ElementResponse(this, 41, Info_pq2);
04197 }
04198
04199
04200 else if (strcmp(argv[0],"GaussPoint") == 0 || strcmp(argv[0],"gausspoint") == 0)
04201 {
04202 return new ElementResponse(this, 6, GaussCoord);
04203 }
04204
04205
04206
04207
04208
04209
04210
04211
04212
04213 else
04214 return 0;
04215 }
04216
04217
04218 int TwentySevenNodeBrick::getResponse (int responseID, Information &eleInfo)
04219 {
04220 switch (responseID) {
04221
04222 case 1:
04223 return eleInfo.setVector( this->getResistingForce() );
04224
04225 case 2:
04226 {
04227
04228 int count = r_integration_order* s_integration_order * t_integration_order;
04229
04230
04231
04232 this->computeGaussPoint();
04233
04234
04235 InfoPlastic(0) = GaussCoord(0);
04236
04237 straintensor pl_stn;
04238
04239 int plastify;
04240 for (int i = 0; i < count; i++) {
04241 plastify = 0;
04242 InfoPlastic(i*4+1) = GaussCoord(i*3+1);
04243 InfoPlastic(i*4+2) = GaussCoord(i*3+2);
04244 InfoPlastic(i*4+3) = GaussCoord(i*3+3);
04245 pl_stn = matpoint[i]->getPlasticStrainTensor();
04246
04247 double q_plastc = pl_stn.q_deviatoric();
04248
04249
04250
04251
04252
04253
04254
04255 InfoPlastic(i*4+4) = q_plastc;
04256
04257 }
04258 return eleInfo.setVector( InfoPlastic );
04259
04260
04261 }
04262 case 3:
04263 {
04264 int count = r_integration_order* s_integration_order * t_integration_order;
04265 stresstensor sts;
04266
04267 this->computeGaussPoint();
04268 Vector wt(9);
04269 int i, rs;
04270
04271
04272 InfoMoment(0) = GaussCoord(0);
04273 Vector Mt(3), Q(3);
04274
04275
04276 InfoMoment(109+0) = GaussCoord(6);
04277
04278
04279 const Vector &coor = theNodes[17]->getCrds();
04280 const Vector &TotDis = theNodes[17]->getTrialDisp();
04281
04282
04283
04284
04285
04286 InfoMoment(109+1) = TotDis(0);
04287
04288
04289
04290
04291 const char *tp = matpoint[1]->getType();
04292 int tag = matpoint[1]->getTag();
04293
04294
04295 float height = 1;
04296
04297 double offset[30];
04298
04299 offset[1] = -0.000; offset[4] = -0.000;
04300 offset[2] = 0.000; offset[5] = 0.000;
04301 offset[3] = 0.000; offset[6] = 0.000;
04302
04303
04304
04305
04306
04307
04308
04309
04310
04311
04312
04313
04314
04315
04316
04317
04318
04319
04320
04321
04322
04323
04324
04325
04326 if (strcmp(tp, "ElasticIsotropic3D") == 0 )
04327 {
04328 wt = getWeightofGP();
04329 const Vector &end1Crd = theNodes[0]->getCrds();
04330 const Vector &end5Crd = theNodes[4]->getCrds();
04331 height = end1Crd(2) - end5Crd(2);
04332
04333
04334
04335
04336 }
04337
04338
04339 Mt(0) = 0; Mt(1) = 0; Mt(2) = 0;
04340
04341
04342
04343 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
04344 {
04345
04346 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
04347 {
04348
04349 rs = (GP_c_r-1)*s_integration_order+GP_c_s-1;
04350
04351 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
04352 {
04353
04354 i =
04355 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
04356
04357 sts = matpoint[i]->getStressTensor();
04358 InfoMoment(i*4+1) = GaussCoord(i*3+1);
04359 InfoMoment(i*4+2) = GaussCoord(i*3+2);
04360 InfoMoment(i*4+3) = GaussCoord(i*3+3);
04361 InfoMoment(i*4+4) = sts.cval(3,3);
04362
04363
04364
04365
04366
04367 if (strcmp(tp, "ElasticIsotropic3D") == 0 ){
04368 Mt(GP_c_t-1) += wt(rs)*sts.cval(3,3)*( InfoMoment(i*4+1)-offset[tag] )/height;
04369
04370
04371 }
04372
04373
04374
04375
04376
04377 }
04378 }
04379 }
04380
04381 InfoMoment(109+2) = ( Mt(0)+Mt(1)+Mt(2) )*0.3333;
04382
04383
04384
04385
04386
04387
04388
04389 return eleInfo.setVector( InfoMoment );
04390 }
04391 case 4:
04392 {
04393 int count = r_integration_order* s_integration_order * t_integration_order;
04394 int i;
04395 stresstensor sts;
04396
04397
04398
04399
04400 InfoStress(0) = count;
04401
04402 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
04403 {
04404
04405 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
04406 {
04407
04408
04409
04410 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
04411 {
04412
04413 i =
04414 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
04415
04416 sts = matpoint[i]->getStressTensor();
04417 InfoStress(i*6+1) = sts.cval(1,1);
04418 InfoStress(i*6+2) = sts.cval(2,2);
04419 InfoStress(i*6+3) = sts.cval(3,3);
04420 InfoStress(i*6+4) = sts.cval(1,2);
04421 InfoStress(i*6+5) = sts.cval(1,3);
04422 InfoStress(i*6+6) = sts.cval(2,3);
04423 }
04424 }
04425 }
04426 return eleInfo.setVector( InfoStress );
04427 }
04428
04429 case 41:
04430 {
04431 int count = r_integration_order* s_integration_order * t_integration_order;
04432 count = count / 2;
04433 stresstensor sts;
04434 sts = matpoint[count]->getStressTensor();
04435 Info_pq2(0) =sts.p_hydrostatic();
04436 Info_pq2(1) =sts.q_deviatoric();
04437 return eleInfo.setVector( Info_pq2 );
04438 }
04439 case 5:
04440 return eleInfo.setMatrix(this->getTangentStiff());
04441
04442 case 6:
04443 {
04444 this->computeGaussPoint();
04445 return eleInfo.setVector( GaussCoord );
04446 }
04447
04448 default:
04449 return -1;
04450 }
04451
04452 }
04453
04454
04455
04457 Vector TwentySevenNodeBrick::getWeightofGP(void)
04458 {
04459
04460
04461
04462
04463 Vector Weight( FixedOrder * FixedOrder );
04464
04465 double r = 0.0;
04466 double rw = 0.0;
04467 double s = 0.0;
04468 double sw = 0.0;
04469 double t = 0.0;
04470 double tw = 0.0;
04471
04472 short where = 0;
04473 short rs = 0;
04474 double tmp = 0;
04475
04476 double weight = 0.0;
04477
04478 int dh_dim[] = {27,3};
04479
04480 tensor dh(2, dh_dim, 0.0);
04481
04482 double det_of_Jacobian = 0.0;
04483
04484 tensor Jacobian;
04485
04486 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
04487 {
04488 r = get_Gauss_p_c( r_integration_order, GP_c_r );
04489 rw = get_Gauss_p_w( r_integration_order, GP_c_r );
04490 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
04491 {
04492 s = get_Gauss_p_c( s_integration_order, GP_c_s );
04493 sw = get_Gauss_p_w( s_integration_order, GP_c_s );
04494
04495 rs = (GP_c_r-1)*s_integration_order+GP_c_s-1;
04496 Weight(rs) = 0;
04497
04498 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
04499 {
04500 t = get_Gauss_p_c( t_integration_order, GP_c_t );
04501 tw = get_Gauss_p_w( t_integration_order, GP_c_t );
04502
04503
04504
04505
04506
04507 dh = dh_drst_at(r,s,t);
04508
04509 Jacobian = Jacobian_3D(dh);
04510 det_of_Jacobian = Jacobian.determinant();
04511
04512
04513
04514 weight = rw * sw * tw * det_of_Jacobian;
04515 Weight(rs) += weight;
04516
04517
04518
04519
04520
04521 }
04522
04523 tmp += Weight(rs);
04524 }
04525 }
04526
04527
04528
04529
04530 return Weight;
04531 }
04532
04533
04535
04536
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546
04547
04548
04549
04550
04551
04552
04553
04554
04555
04556
04557
04558
04559
04560
04561
04562
04563
04564
04565
04566
04567
04568
04569
04570
04571
04572
04573
04574
04575
04576
04577
04578
04579
04580
04581
04582
04583
04584
04585
04586
04587
04588
04589
04590
04591
04592
04593
04594
04595
04596
04597
04598
04599
04600
04601
04602
04603
04604
04605
04606
04607
04608
04609
04610
04611
04612
04613
04614
04615
04616
04617
04618
04619
04620
04621
04622
04623
04624
04625
04626
04627
04628
04629
04630
04631
04632
04633
04634
04635
04636
04637
04638
04639
04640
04641
04642
04643
04644
04645
04646
04647
04648
04649
04650
04651
04652
04653
04654
04655
04656
04657
04658
04659
04660
04661
04662
04663
04664
04665
04666
04667
04668
04669
04670
04671
04672
04673
04674
04675
04676
04677
04678
04679
04680
04681
04682
04683
04684
04685
04686
04687
04688
04689
04690
04691
04692
04693
04694
04695
04696
04697
04698
04699
04700
04701
04702
04703
04704
04705
04706
04707
04708
04709
04710
04711
04712
04713
04714
04715
04716
04717
04718
04719
04720
04721
04722
04723
04724
04725
04726
04727
04728
04729
04730
04731
04732
04733
04734
04735
04736
04737
04738
04739
04740
04741
04742
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781
04782
04783
04784
04785
04786
04787
04788
04789
04790
04791
04792
04793
04794
04795
04796
04798
04799
04800
04801
04802
04803
04804
04805
04806
04807
04808
04809
04810
04811
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822
04823
04824
04825
04826
04827
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838
04839
04840
04841
04842
04843
04844
04845
04846
04847
04848
04849
04851
04852
04853
04858
04860
04861
04862
04863
04864
04865
04866
04867
04868
04869
04870
04871
04872
04874
04875
04876
04877
04878
04879
04880
04881
04882
04883
04884
04888
04889
04890
04891
04892
04893
04894
04895
04896
04897
04898
04899
04900
04901
04902
04903
04904
04905
04906
04907
04908
04909
04910
04911
04912
04913
04914
04915
04916
04917
04918
04919
04920
04921
04922
04923
04924
04925
04926
04927
04928
04929
04930
04931
04932
04933
04934
04935
04936
04937
04938
04939
04940
04941
04942
04943
04944
04945
04946
04947
04948
04949
04950
04951
04952
04953
04954
04955
04956
04957
04958
04959
04960
04961
04962
04963
04964
04965
04966
04967
04968
04969
04970
04971
04973
04974
04975
04976
04977
04978
04979
04980
04981
04982
04983
04984
04985
04986
04987
04988
04989
04990
04991
04992
04993
04994
04995
04996
04997
04998
04999
05000
05001
05002
05003
05004
05005
05006
05007
05008
05009
05012
05013
05014
05015
05016
05017
05018
05019
05020
05021
05022
05023
05024
05025
05026
05027
05028
05029
05030
05031
05032
05033
05034
05035
05036
05037
05038
05039
05040
05041
05042
05043
05045
05046
05047
05048
05049
05050
05051
05052
05053
05054
05055
05056
05057
05058
05059
05060
05061
05062
05063
05064
05065
05066
05067
05068
05069
05070
05071
05072
05073
05074
05075
05076
05077
05078
05079
05080
05081
05082
05083
05084
05085
05086
05087
05088
05089
05090
05091
05092
05093
05094
05095
05096
05097
05098
05099
05100
05101
05102
05103
05104
05105
05106
05107
05108
05109
05110
05111
05112
05113
05114
05115
05116
05117
05118
05119
05120
05121
05122
05123
05124
05125
05126
05127
05128
05129
05130
05131
05132
05133
05134
05135
05136
05137
05138 double TwentySevenNodeBrick::get_Gauss_p_c(short order, short point_numb)
05139 {
05140
05141
05142 static double Gauss_coordinates[7][7];
05143
05144 Gauss_coordinates[1][1] = 0.0 ;
05145 Gauss_coordinates[2][1] = -0.577350269189626;
05146 Gauss_coordinates[2][2] = -Gauss_coordinates[2][1];
05147 Gauss_coordinates[3][1] = -0.774596669241483;
05148 Gauss_coordinates[3][2] = 0.0;
05149 Gauss_coordinates[3][3] = -Gauss_coordinates[3][1];
05150 Gauss_coordinates[4][1] = -0.861136311594053;
05151 Gauss_coordinates[4][2] = -0.339981043584856;
05152 Gauss_coordinates[4][3] = -Gauss_coordinates[4][2];
05153 Gauss_coordinates[4][4] = -Gauss_coordinates[4][1];
05154 Gauss_coordinates[5][1] = -0.906179845938664;
05155 Gauss_coordinates[5][2] = -0.538469310105683;
05156 Gauss_coordinates[5][3] = 0.0;
05157 Gauss_coordinates[5][4] = -Gauss_coordinates[5][2];
05158 Gauss_coordinates[5][5] = -Gauss_coordinates[5][1];
05159 Gauss_coordinates[6][1] = -0.932469514203152;
05160 Gauss_coordinates[6][2] = -0.661509386466265;
05161 Gauss_coordinates[6][3] = -0.238619188183197;
05162 Gauss_coordinates[6][4] = -Gauss_coordinates[6][3];
05163 Gauss_coordinates[6][5] = -Gauss_coordinates[6][2];
05164 Gauss_coordinates[6][6] = -Gauss_coordinates[6][1];
05165
05166 return Gauss_coordinates[order][point_numb];
05167 }
05168
05169 double TwentySevenNodeBrick::get_Gauss_p_w(short order, short point_numb)
05170 {
05171
05172
05173 static double Gauss_weights[7][7];
05174
05175 Gauss_weights[1][1] = 2.0;
05176 Gauss_weights[2][1] = 1.0;
05177 Gauss_weights[2][2] = 1.0;
05178 Gauss_weights[3][1] = 0.555555555555556;
05179 Gauss_weights[3][2] = 0.888888888888889;
05180 Gauss_weights[3][3] = Gauss_weights[3][1];
05181 Gauss_weights[4][1] = 0.347854845137454;
05182 Gauss_weights[4][2] = 0.652145154862546;
05183 Gauss_weights[4][3] = Gauss_weights[4][2];
05184 Gauss_weights[4][4] = Gauss_weights[4][1];
05185 Gauss_weights[5][1] = 0.236926885056189;
05186 Gauss_weights[5][2] = 0.478628670499366;
05187 Gauss_weights[5][3] = 0.568888888888889;
05188 Gauss_weights[5][4] = Gauss_weights[5][2];
05189 Gauss_weights[5][5] = Gauss_weights[5][1];
05190 Gauss_weights[6][1] = 0.171324492379170;
05191 Gauss_weights[6][2] = 0.381761573048139;
05192 Gauss_weights[6][3] = 0.467913934572691;
05193 Gauss_weights[6][4] = Gauss_weights[6][3];
05194 Gauss_weights[6][5] = Gauss_weights[6][2];
05195 Gauss_weights[6][6] = Gauss_weights[6][1];
05196
05197 return Gauss_weights[order][point_numb];
05198 }
05199
05200 int TwentySevenNodeBrick::update()
05201 {
05202 double r = 0.0;
05203
05204 double s = 0.0;
05205
05206 double t = 0.0;
05207
05208
05209 short where = 0;
05210
05211
05212 int dh_dim[] = {27,3};
05213 tensor dh(2, dh_dim, 0.0);
05214
05215
05216 static int disp_dim[] = {27,3};
05217 tensor incremental_displacements(2,disp_dim,0.0);
05218
05219 straintensor incremental_strain;
05220
05221 tensor Jacobian;
05222 tensor JacobianINV;
05223 tensor dhGlobal;
05224
05225 incremental_displacements = incr_disp();
05226
05227 for( short GP_c_r = 1 ; GP_c_r <= r_integration_order ; GP_c_r++ )
05228 {
05229 r = get_Gauss_p_c( r_integration_order, GP_c_r );
05230
05231 for( short GP_c_s = 1 ; GP_c_s <= s_integration_order ; GP_c_s++ )
05232 {
05233 s = get_Gauss_p_c( s_integration_order, GP_c_s );
05234
05235 for( short GP_c_t = 1 ; GP_c_t <= t_integration_order ; GP_c_t++ )
05236 {
05237 t = get_Gauss_p_c( t_integration_order, GP_c_t );
05238
05239
05240
05241 where =
05242 ((GP_c_r-1)*s_integration_order+GP_c_s-1)*t_integration_order+GP_c_t-1;
05243
05244 dh = dh_drst_at(r,s,t);
05245
05246 Jacobian = Jacobian_3D(dh);
05247
05248
05249 JacobianINV = Jacobian_3Dinv(dh);
05250
05251
05252
05253
05254
05255
05256 dhGlobal = dh("ij") * JacobianINV("kj");
05257
05258
05259
05260
05261
05262
05263
05264
05265
05266
05267
05268
05269
05270 incremental_strain =
05271 (dhGlobal("ib")*incremental_displacements("ia")).symmetrize11();
05272 incremental_strain.null_indices();
05273
05274
05275
05276
05277
05278
05279 if ( ( (matpoint[where]->matmodel)->setTrialStrainIncr( incremental_strain)) )
05280 opserr << "TwentySevenNodeBrick::update (tag: " << this->getTag() << "), update() failed\n";
05281 }
05282 }
05283 }
05284 return 0;
05285 }
05286
05287
05288
05289 #endif