00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <math.h>
00034
00035 #include <ID.h>
00036 #include <Vector.h>
00037 #include <Matrix.h>
00038 #include <Element.h>
00039 #include <Node.h>
00040 #include <NDMaterial.h>
00041 #include <Domain.h>
00042 #include <ErrorHandler.h>
00043 #include <NineNodeMixedQuad.h>
00044 #include <ElementResponse.h>
00045
00046 #include <Renderer.h>
00047
00048 #include <Channel.h>
00049 #include <FEM_ObjectBroker.h>
00050
00051
00052 Matrix NineNodeMixedQuad :: stiff(18,18) ;
00053 Vector NineNodeMixedQuad :: resid(18) ;
00054 Matrix NineNodeMixedQuad :: mass(18,18) ;
00055 double NineNodeMixedQuad::xl[2][9];
00056
00057
00058 double NineNodeMixedQuad::root06 = sqrt(0.6) ;
00059 double NineNodeMixedQuad::sg[] = { -root06, 0.0, root06 } ;
00060 double NineNodeMixedQuad::wg[] = { 5.0/9.0, 8.0/9.0, 5.0/9.0 } ;
00061
00062
00063
00064 NineNodeMixedQuad :: NineNodeMixedQuad( ) :
00065 Element( 0, ELE_TAG_NineNodeMixedQuad ),
00066 connectedExternalNodes(9) , load(0), Ki(0)
00067 {
00068 for (int i=0; i<9; i++)
00069 materialPointers[i] = 0;
00070 }
00071
00072
00073
00074 NineNodeMixedQuad :: NineNodeMixedQuad( int tag,
00075 int node1,
00076 int node2,
00077 int node3,
00078 int node4,
00079 int node5,
00080 int node6,
00081 int node7,
00082 int node8,
00083 int node9,
00084 NDMaterial &theMaterial ) :
00085 Element( tag, ELE_TAG_NineNodeMixedQuad ),
00086 connectedExternalNodes(9) , load(0), Ki(0)
00087 {
00088 connectedExternalNodes(0) = node1 ;
00089 connectedExternalNodes(1) = node2 ;
00090 connectedExternalNodes(2) = node3 ;
00091 connectedExternalNodes(3) = node4 ;
00092 connectedExternalNodes(4) = node5 ;
00093 connectedExternalNodes(5) = node6 ;
00094 connectedExternalNodes(6) = node7 ;
00095 connectedExternalNodes(7) = node8 ;
00096 connectedExternalNodes(8) = node9 ;
00097
00098 int i ;
00099 for ( i=0 ; i<9; i++ ) {
00100
00101 materialPointers[i] = theMaterial.getCopy("AxiSymmetric2D") ;
00102
00103 if (materialPointers[i] == 0) {
00104
00105 opserr << "NineNodeMixedQuad::constructor() - failed to get a material of type: AxiSymmetric2D\n";
00106 }
00107
00108 }
00109 }
00110
00111
00112
00113 NineNodeMixedQuad :: ~NineNodeMixedQuad( )
00114 {
00115 int i ;
00116 for ( i=0 ; i<9; i++ ) {
00117
00118 delete materialPointers[i] ;
00119 materialPointers[i] = 0 ;
00120
00121 nodePointers[i] = 0 ;
00122 }
00123
00124 if (load != 0)
00125 delete load;
00126
00127 if (Ki != 0)
00128 delete Ki;
00129 }
00130
00131
00132
00133 void
00134 NineNodeMixedQuad::setDomain( Domain *theDomain )
00135 {
00136 for ( int i = 0; i<9; i++ )
00137 nodePointers[i] = theDomain->getNode( connectedExternalNodes(i) ) ;
00138
00139 this->DomainComponent::setDomain(theDomain);
00140 }
00141
00142
00143
00144 int
00145 NineNodeMixedQuad::getNumExternalNodes( ) const
00146 {
00147 return 9 ;
00148 }
00149
00150
00151
00152 const ID&
00153 NineNodeMixedQuad::getExternalNodes( )
00154 {
00155 return connectedExternalNodes ;
00156 }
00157
00158 Node **
00159 NineNodeMixedQuad::getNodePtrs(void)
00160 {
00161 return nodePointers;
00162 }
00163
00164
00165
00166 int
00167 NineNodeMixedQuad::getNumDOF( )
00168 {
00169 return 18 ;
00170 }
00171
00172
00173
00174 int
00175 NineNodeMixedQuad::commitState( )
00176 {
00177 int i ;
00178 int success = 0 ;
00179
00180
00181 if ((success = this->Element::commitState()) != 0) {
00182 opserr << "NineNodeMixedQuad::commitState () - failed in base class\n";
00183 }
00184
00185 for ( i=0; i<9; i++ )
00186 success += materialPointers[i]->commitState( ) ;
00187
00188 return success ;
00189 }
00190
00191
00192
00193
00194 int
00195 NineNodeMixedQuad::revertToLastCommit( )
00196 {
00197 int i ;
00198 int success = 0 ;
00199
00200 for ( i=0; i<9; i++ )
00201 success += materialPointers[i]->revertToLastCommit( ) ;
00202
00203 return success ;
00204 }
00205
00206
00207
00208 int
00209 NineNodeMixedQuad::revertToStart( )
00210 {
00211 int i ;
00212 int success = 0 ;
00213
00214 for ( i=0; i<9; i++ )
00215 success += materialPointers[i]->revertToStart( ) ;
00216
00217 return success ;
00218 }
00219
00220
00221 void
00222 NineNodeMixedQuad::Print( OPS_Stream &s, int flag )
00223 {
00224 s << endln ;
00225 s << "Nine Node Quad -- Mixed Pressure/Volume -- Plane Strain \n" ;
00226 s << "Element Number " << this->getTag() << endln ;
00227 s << "Node 1 : " << connectedExternalNodes(0) << endln ;
00228 s << "Node 2 : " << connectedExternalNodes(1) << endln ;
00229 s << "Node 3 : " << connectedExternalNodes(2) << endln ;
00230 s << "Node 4 : " << connectedExternalNodes(3) << endln ;
00231 s << "Node 5 : " << connectedExternalNodes(4) << endln ;
00232 s << "Node 6 : " << connectedExternalNodes(5) << endln ;
00233 s << "Node 7 : " << connectedExternalNodes(6) << endln ;
00234 s << "Node 8 : " << connectedExternalNodes(7) << endln ;
00235 s << "Node 9 : " << connectedExternalNodes(8) << endln ;
00236
00237 s << "Material Information : \n " ;
00238
00239 materialPointers[0]->Print( s, flag ) ;
00240
00241 s << endln ;
00242 }
00243
00244
00245
00246 const Matrix&
00247 NineNodeMixedQuad::getTangentStiff( )
00248 {
00249 int tang_flag = 1 ;
00250
00251
00252 formResidAndTangent( tang_flag ) ;
00253
00254 return stiff ;
00255 }
00256
00257
00258
00259 const Matrix&
00260 NineNodeMixedQuad::getInitialStiff( )
00261 {
00262 if (Ki != 0)
00263 return *Ki;
00264
00265 static const int ndm = 2 ;
00266
00267 static const int ndf = 2 ;
00268
00269 static const int nstress = 4 ;
00270
00271 static const int numberNodes = 9 ;
00272
00273 static const int numberGauss = 9 ;
00274
00275 static const int nShape = 3 ;
00276
00277 static const int nMixed = 3 ;
00278
00279 int i, j, k, p, q, r, s ;
00280 int jj, kk ;
00281
00282 static double volume ;
00283
00284 static double xsj ;
00285
00286 static double dvol[numberGauss] ;
00287
00288 static double gaussPoint[ndm] ;
00289
00290 static double natCoorArray[ndm][numberGauss] ;
00291
00292 static Vector strain(nstress) ;
00293
00294 static double shp[nShape][numberNodes] ;
00295
00296 static double Shape[nShape][numberNodes][numberGauss] ;
00297
00298 static double shpBar[nShape][numberNodes][nMixed] ;
00299
00300 static double rightHandSide[nShape][numberNodes][nMixed] ;
00301
00302 static Vector residJ(ndf) ;
00303
00304 static Matrix stiffJK(ndf,ndf) ;
00305
00306 static Vector stress(nstress) ;
00307
00308 static Matrix dd(nstress,nstress) ;
00309
00310 static double interp[nMixed] ;
00311
00312 static Matrix Proj(3,3) ;
00313 static Matrix ProjInv(3,3) ;
00314
00315 static Matrix Iden(3,3) ;
00316 Iden(0,0) = 1.0 ;
00317 Iden(1,1) = 1.0 ;
00318 Iden(2,2) = 1.0 ;
00319
00320
00321
00322 static Matrix BJ(nstress,ndf) ;
00323
00324 static Matrix BJtran(ndf,nstress) ;
00325
00326 static Matrix BK(nstress,ndf) ;
00327
00328 static Matrix BJtranD(ndf,nstress) ;
00329
00330
00331
00332
00333
00334 stiff.Zero( ) ;
00335
00336
00337 computeBasis() ;
00338
00339
00340 for ( p=0; p<nShape; p++ ) {
00341 for ( q=0; q<numberNodes; q++ ) {
00342
00343 for (r=0; r<nMixed; r++ ) {
00344 shpBar[p][q][r] = 0.0 ;
00345 rightHandSide[p][q][r] = 0.0 ;
00346 }
00347
00348 }
00349 }
00350
00351
00352
00353 volume = 0.0 ;
00354
00355
00356 Proj.Zero( ) ;
00357 ProjInv.Zero( ) ;
00358
00359
00360 int count = 0 ;
00361
00362 for ( i = 0; i < 3; i++ ) {
00363 for ( j = 0; j < 3; j++ ) {
00364
00365 gaussPoint[0] = sg[i] ;
00366 gaussPoint[1] = sg[j] ;
00367
00368
00369
00370 natCoorArray[0][count] = gaussPoint[0] ;
00371 natCoorArray[1][count] = gaussPoint[1] ;
00372
00373
00374
00375 shape2dNine( gaussPoint, xl, shp, xsj ) ;
00376
00377
00378
00379 for ( p=0; p<nShape; p++ ) {
00380 for ( q=0; q<numberNodes; q++ )
00381 Shape[p][q][count] = shp[p][q] ;
00382 }
00383
00384
00385
00386 dvol[count] = ( wg[i]*wg[j] ) * xsj ;
00387
00388
00389
00390 interp[0] = 1.0 ;
00391 interp[1] = gaussPoint[0] ;
00392 interp[2] = gaussPoint[1] ;
00393
00394 for ( r=0; r<nMixed; r++ ) {
00395 for ( s=0; s<nMixed; s++ )
00396 Proj(r,s) += ( interp[r]*interp[s] * dvol[count] ) ;
00397 }
00398
00399 volume += dvol[count] ;
00400
00401
00402
00403 for ( p=0; p<nShape; p++ ) {
00404 for ( q=0; q<numberNodes; q++ ) {
00405
00406 for ( s=0; s<nMixed; s++ )
00407 rightHandSide[p][q][s] += ( shp[p][q] * interp[s] * dvol[count] ) ;
00408
00409 }
00410 }
00411
00412
00413
00414 count++ ;
00415
00416 }
00417 }
00418
00419
00420
00421
00422
00423 Proj.Solve( Iden, ProjInv ) ;
00424
00425
00426 for ( p=0; p<nShape; p++ ) {
00427 for ( q=0; q<numberNodes; q++ ) {
00428
00429 for (r=0; r<nMixed; r++ ) {
00430 for (s=0; s<nMixed; s++ )
00431 shpBar[p][q][r] += ( ProjInv(r,s) * rightHandSide[p][q][s] ) ;
00432 }
00433
00434 }
00435 }
00436
00437
00438
00439 for ( i=0; i<numberGauss; i++ ) {
00440
00441
00442 gaussPoint[0] = natCoorArray[0][i] ;
00443 gaussPoint[1] = natCoorArray[1][i] ;
00444
00445
00446 for ( p=0; p<nShape; p++ ) {
00447 for ( q=0; q<numberNodes; q++ )
00448 shp[p][q] = Shape[p][q][i] ;
00449 }
00450
00451 dd = materialPointers[i]->getInitialTangent( ) ;
00452 dd *= dvol[i] ;
00453
00454
00455
00456
00457 jj = 0 ;
00458 for ( j=0; j<numberNodes; j++ ) {
00459
00460 BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00461
00462
00463
00464 for (p=0; p<ndf; p++) {
00465 for (q=0; q<nstress; q++)
00466 BJtran(p,q) = BJ(q,p) ;
00467 }
00468
00469
00470
00471 BJtranD.addMatrixProduct(0.0, BJtran,dd,1.0);
00472
00473 kk = 0 ;
00474 for ( k=0; k<numberNodes; k++ ) {
00475
00476 BK = computeBbar( k, gaussPoint, shp, shpBar ) ;
00477
00478
00479
00480 stiffJK.addMatrixProduct(0.0, BJtranD,BK,1.0) ;
00481
00482 for ( p=0; p<ndf; p++ ) {
00483 for ( q=0; q<ndf; q++ )
00484 stiff( jj+p, kk+q ) += stiffJK( p, q ) ;
00485 }
00486
00487 kk += ndf ;
00488 }
00489
00490 jj += ndf ;
00491 }
00492 }
00493
00494 Ki = new Matrix(stiff);
00495
00496 return stiff;
00497 }
00498
00499
00500
00501 const Matrix&
00502 NineNodeMixedQuad::getMass( )
00503 {
00504 int tangFlag = 1 ;
00505
00506 formInertiaTerms( tangFlag ) ;
00507
00508 return mass ;
00509 }
00510
00511
00512
00513 void
00514 NineNodeMixedQuad::zeroLoad( )
00515 {
00516 if (load != 0)
00517 load->Zero();
00518 }
00519
00520
00521
00522 int
00523 NineNodeMixedQuad::addLoad(ElementalLoad *theLoad, double loadFactor)
00524 {
00525 opserr << "NineNodeMixedQuad::addLoad - load type unknown for ele with tag: " << this->getTag() << endln;
00526 return -1;
00527 }
00528
00529 int
00530 NineNodeMixedQuad::addInertiaLoadToUnbalance(const Vector &accel)
00531 {
00532 static const int numberGauss = 9 ;
00533 static const int numberNodes = 9 ;
00534 static const int ndf = 2 ;
00535
00536 int i;
00537
00538
00539 int haveRho = 0;
00540 for (i = 0; i < numberGauss; i++) {
00541 if (materialPointers[i]->getRho() != 0.0)
00542 haveRho = 1;
00543 }
00544
00545 if (haveRho == 0)
00546 return 0;
00547
00548
00549 int tangFlag = 1 ;
00550 formInertiaTerms( tangFlag ) ;
00551
00552
00553 int count = 0;
00554
00555 for (i=0; i<numberNodes; i++) {
00556 const Vector &Raccel = nodePointers[i]->getRV(accel);
00557 for (int j=0; j<ndf; j++)
00558 resid(count++) = Raccel(i);
00559 }
00560
00561
00562 if (load == 0)
00563 load = new Vector(numberNodes*ndf);
00564
00565
00566 load->addMatrixVector(1.0, mass, resid, -1.0);
00567
00568 return 0;
00569 }
00570
00571
00572
00573 const Vector&
00574 NineNodeMixedQuad::getResistingForce( )
00575 {
00576 int tang_flag = 0 ;
00577
00578 formResidAndTangent( tang_flag ) ;
00579
00580
00581 if (load != 0)
00582 resid -= *load;
00583
00584 return resid ;
00585 }
00586
00587
00588
00589 const Vector&
00590 NineNodeMixedQuad::getResistingForceIncInertia( )
00591 {
00592 int tang_flag = 0 ;
00593
00594 static Vector res(18);
00595
00596
00597 formResidAndTangent( tang_flag ) ;
00598
00599
00600 formInertiaTerms( tang_flag ) ;
00601
00602 res = resid;
00603
00604
00605 if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
00606 res += this->getRayleighDampingForces();
00607
00608
00609 if (load != 0)
00610 res -= *load;
00611
00612 return res;
00613 }
00614
00615
00616
00617
00618 void
00619 NineNodeMixedQuad::formInertiaTerms( int tangFlag )
00620 {
00621
00622 static const int ndm = 2 ;
00623
00624 static const int ndf = 2 ;
00625
00626 static const int numberNodes = 9 ;
00627
00628
00629
00630 static const int nShape = 3 ;
00631
00632 static const int massIndex = nShape - 1 ;
00633
00634 double xsj ;
00635
00636 double dvol ;
00637
00638 static double shp[nShape][numberNodes] ;
00639
00640 static Vector momentum(ndf) ;
00641
00642 static Matrix sx(ndm,ndm) ;
00643
00644 static double GaussPoint[2] ;
00645
00646 int j, k, p, q, r ;
00647 int jj, kk ;
00648
00649 double temp, rho, massJK ;
00650
00651
00652
00653 mass.Zero( ) ;
00654
00655
00656 computeBasis() ;
00657
00658
00659 int count = 0 ;
00660 for ( p=0; p<3; p++ ) {
00661 for ( q=0; q<3; q++ ) {
00662
00663 GaussPoint[0] = sg[p] ;
00664 GaussPoint[1] = sg[q] ;
00665
00666
00667 shape2dNine( GaussPoint, xl, shp, xsj ) ;
00668
00669
00670 dvol = ( wg[p] * wg[q] ) * xsj ;
00671
00672
00673 momentum.Zero( ) ;
00674 for ( j = 0; j < numberNodes; j++ )
00675
00676 momentum.addVector( 1.0,
00677 nodePointers[j]->getTrialAccel(),
00678 shp[massIndex][j] ) ;
00679
00680
00681
00682 rho = materialPointers[count]->getRho() ;
00683
00684
00685 momentum *= rho ;
00686
00687
00688
00689 for ( jj=0, j=0; j<numberNodes; j++, jj+=ndf ) {
00690
00691 temp = shp[massIndex][j] * dvol ;
00692
00693 for ( r=0; r<ndf; r++ )
00694 resid( jj+r ) += ( temp * momentum(r) ) ;
00695
00696
00697 if ( tangFlag == 1 ) {
00698
00699
00700 temp *= rho ;
00701
00702
00703 for ( kk=0, k=0; k<numberNodes; k++, kk+=ndf ) {
00704
00705 massJK = temp * shp[massIndex][k] ;
00706
00707 for ( r=0; r<ndf; r++ )
00708 mass( jj+r, kk+r ) += massJK ;
00709
00710 }
00711
00712 }
00713
00714 }
00715
00716
00717 count++ ;
00718 }
00719 }
00720
00721 }
00722
00723
00724
00725
00726 void
00727 NineNodeMixedQuad::formResidAndTangent( int tang_flag )
00728 {
00729
00730
00731
00732
00733 static const int ndm = 2 ;
00734
00735 static const int ndf = 2 ;
00736
00737 static const int nstress = 4 ;
00738
00739 static const int numberNodes = 9 ;
00740
00741 static const int numberGauss = 9 ;
00742
00743 static const int nShape = 3 ;
00744
00745 static const int nMixed = 3 ;
00746
00747 int i, j, k, p, q, r, s ;
00748 int jj, kk ;
00749
00750 int success ;
00751
00752 static double volume ;
00753
00754 static double xsj ;
00755
00756 static double dvol[numberGauss] ;
00757
00758 static double gaussPoint[ndm] ;
00759
00760 static double natCoorArray[ndm][numberGauss] ;
00761
00762 static Vector strain(nstress) ;
00763
00764 static double shp[nShape][numberNodes] ;
00765
00766 static double Shape[nShape][numberNodes][numberGauss] ;
00767
00768 static double shpBar[nShape][numberNodes][nMixed] ;
00769
00770 static double rightHandSide[nShape][numberNodes][nMixed] ;
00771
00772 static Vector residJ(ndf) ;
00773
00774 static Matrix stiffJK(ndf,ndf) ;
00775
00776 static Vector stress(nstress) ;
00777
00778 static Matrix dd(nstress,nstress) ;
00779
00780 static double interp[nMixed] ;
00781
00782 static Matrix Proj(3,3) ;
00783 static Matrix ProjInv(3,3) ;
00784
00785 static Matrix Iden(3,3) ;
00786 Iden(0,0) = 1.0 ;
00787 Iden(1,1) = 1.0 ;
00788 Iden(2,2) = 1.0 ;
00789
00790
00791
00792 static Matrix BJ(nstress,ndf) ;
00793
00794 static Matrix BJtran(ndf,nstress) ;
00795
00796 static Matrix BK(nstress,ndf) ;
00797
00798 static Matrix BJtranD(ndf,nstress) ;
00799
00800
00801
00802
00803
00804 stiff.Zero( ) ;
00805 resid.Zero( ) ;
00806
00807
00808 computeBasis() ;
00809
00810
00811 for ( p=0; p<nShape; p++ ) {
00812 for ( q=0; q<numberNodes; q++ ) {
00813
00814 for (r=0; r<nMixed; r++ ) {
00815 shpBar[p][q][r] = 0.0 ;
00816 rightHandSide[p][q][r] = 0.0 ;
00817 }
00818
00819 }
00820 }
00821
00822
00823
00824 volume = 0.0 ;
00825
00826
00827 Proj.Zero( ) ;
00828 ProjInv.Zero( ) ;
00829
00830
00831
00832 int count = 0 ;
00833
00834 for ( i = 0; i < 3; i++ ) {
00835 for ( j = 0; j < 3; j++ ) {
00836
00837 gaussPoint[0] = sg[i] ;
00838 gaussPoint[1] = sg[j] ;
00839
00840
00841
00842 natCoorArray[0][count] = gaussPoint[0] ;
00843 natCoorArray[1][count] = gaussPoint[1] ;
00844
00845
00846
00847 shape2dNine( gaussPoint, xl, shp, xsj ) ;
00848
00849
00850
00851 for ( p=0; p<nShape; p++ ) {
00852 for ( q=0; q<numberNodes; q++ )
00853 Shape[p][q][count] = shp[p][q] ;
00854 }
00855
00856
00857
00858 dvol[count] = ( wg[i]*wg[j] ) * xsj ;
00859
00860
00861
00862 interp[0] = 1.0 ;
00863 interp[1] = gaussPoint[0] ;
00864 interp[2] = gaussPoint[1] ;
00865
00866 for ( r=0; r<nMixed; r++ ) {
00867 for ( s=0; s<nMixed; s++ )
00868 Proj(r,s) += ( interp[r]*interp[s] * dvol[count] ) ;
00869 }
00870
00871 volume += dvol[count] ;
00872
00873
00874
00875 for ( p=0; p<nShape; p++ ) {
00876 for ( q=0; q<numberNodes; q++ ) {
00877
00878 for ( s=0; s<nMixed; s++ )
00879 rightHandSide[p][q][s] += ( shp[p][q] * interp[s] * dvol[count] ) ;
00880
00881 }
00882 }
00883
00884
00885
00886 count++ ;
00887
00888 }
00889 }
00890
00891
00892
00893
00894
00895 Proj.Solve( Iden, ProjInv ) ;
00896
00897
00898 for ( p=0; p<nShape; p++ ) {
00899 for ( q=0; q<numberNodes; q++ ) {
00900
00901 for (r=0; r<nMixed; r++ ) {
00902 for (s=0; s<nMixed; s++ )
00903 shpBar[p][q][r] += ( ProjInv(r,s) * rightHandSide[p][q][s] ) ;
00904 }
00905
00906 }
00907 }
00908
00909
00910
00911 for ( i=0; i<numberGauss; i++ ) {
00912
00913
00914 gaussPoint[0] = natCoorArray[0][i] ;
00915 gaussPoint[1] = natCoorArray[1][i] ;
00916
00917
00918 for ( p=0; p<nShape; p++ ) {
00919 for ( q=0; q<numberNodes; q++ )
00920 shp[p][q] = Shape[p][q][i] ;
00921 }
00922
00923
00924
00925 strain.Zero( ) ;
00926
00927
00928 for ( j=0; j<numberNodes; j++ ) {
00929
00930
00931
00932 BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00933
00934
00935 const Vector &ul = nodePointers[j]->getTrialDisp( ) ;
00936
00937
00938
00939 strain.addMatrixVector(1.0, BJ,ul,1.0 ) ;
00940
00941 }
00942
00943
00944
00945
00946 success = materialPointers[i]->setTrialStrain( strain ) ;
00947
00948
00949 stress = materialPointers[i]->getStress( ) ;
00950
00951
00952
00953 stress *= dvol[i] ;
00954
00955 if ( tang_flag == 1 ) {
00956 dd = materialPointers[i]->getTangent( ) ;
00957 dd *= dvol[i] ;
00958 }
00959
00960
00961
00962
00963 jj = 0 ;
00964 for ( j=0; j<numberNodes; j++ ) {
00965
00966 BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00967
00968
00969
00970 for (p=0; p<ndf; p++) {
00971 for (q=0; q<nstress; q++)
00972 BJtran(p,q) = BJ(q,p) ;
00973 }
00974
00975
00976
00977
00978 residJ.addMatrixVector(0.0, BJtran,stress,1.0);
00979
00980
00981 for ( p=0; p<ndf; p++ )
00982 resid( jj + p ) += residJ(p) ;
00983
00984
00985 if ( tang_flag == 1 ) {
00986
00987
00988 BJtranD.addMatrixProduct(0.0, BJtran,dd,1.0);
00989
00990 kk = 0 ;
00991 for ( k=0; k<numberNodes; k++ ) {
00992
00993 BK = computeBbar( k, gaussPoint, shp, shpBar ) ;
00994
00995
00996
00997 stiffJK.addMatrixProduct(0.0, BJtranD,BK,1.0) ;
00998
00999 for ( p=0; p<ndf; p++ ) {
01000 for ( q=0; q<ndf; q++ )
01001 stiff( jj+p, kk+q ) += stiffJK( p, q ) ;
01002 }
01003
01004 kk += ndf ;
01005 }
01006
01007 }
01008
01009 jj += ndf ;
01010 }
01011
01012
01013 }
01014
01015
01016 return ;
01017 }
01018
01019
01020
01021
01022
01023 void
01024 NineNodeMixedQuad::computeBasis( )
01025 {
01026
01027
01028
01029 for ( int i = 0; i<9; i++ ) {
01030
01031 const Vector &coorI = nodePointers[i]->getCrds( ) ;
01032
01033 xl[0][i] = coorI(0) ;
01034 xl[1][i] = coorI(1) ;
01035
01036 }
01037
01038 }
01039
01040
01041
01042
01043 const Matrix&
01044 NineNodeMixedQuad::computeBbar( int node,
01045 const double natCoor[2],
01046 const double shp[3][9],
01047 double shpBar[3][9][3] )
01048 {
01049
01050 static Matrix Bbar(4,2) ;
01051
01052 static double Bdev[3][2] ;
01053
01054 static double BbarVol[3][2] ;
01055
01056 static const double one3 = 1.0/3.0 ;
01057
01058 static double interp[3] ;
01059
01060 static double c0, c1 ;
01061
01062 int i, j ;
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 Bbar.Zero( ) ;
01090
01091
01092 Bdev[0][0] = 2.0*shp[0][node] ;
01093 Bdev[0][1] = -shp[1][node] ;
01094
01095 Bdev[1][0] = -shp[0][node] ;
01096 Bdev[1][1] = 2.0*shp[1][node] ;
01097
01098 Bdev[2][0] = -shp[0][node] ;
01099 Bdev[2][1] = -shp[1][node] ;
01100
01101
01102
01103 interp[0] = 1.0 ;
01104 interp[1] = natCoor[0] ;
01105 interp[2] = natCoor[1] ;
01106
01107
01108 c0 = 0.0 ;
01109 c1 = 0.0 ;
01110
01111 for (i=0; i<3; i++) {
01112
01113 c0 += ( shpBar[0][node][i] * interp[i] ) ;
01114 c1 += ( shpBar[1][node][i] * interp[i] ) ;
01115
01116 }
01117
01118
01119
01120
01121
01122 BbarVol[0][0] = c0 ;
01123 BbarVol[0][1] = c1 ;
01124
01125 BbarVol[1][0] = c0 ;
01126 BbarVol[1][1] = c1 ;
01127
01128 BbarVol[2][0] = c0 ;
01129 BbarVol[2][1] = c1 ;
01130
01131
01132
01133 for ( i=0; i<3; i++ ){
01134 for ( j=0; j<2; j++ )
01135 Bbar(i,j) = one3*( Bdev[i][j] + BbarVol[i][j] ) ;
01136 }
01137
01138
01139
01140 Bbar(3,0) = shp[1][node] ;
01141 Bbar(3,1) = shp[0][node] ;
01142
01143
01144 return Bbar ;
01145 }
01146
01147
01148
01149
01150 void
01151 NineNodeMixedQuad::shape2dNine( double coor[2],
01152 const double x[2][9],
01153 double shp[3][9],
01154 double &xsj )
01155 {
01156 static const int nNode = 9 ;
01157 static const int ndm = 2 ;
01158
01159 static const int node1[] = { 0,1,1,0, 2,1,2,0, 2 } ;
01160 static const int node2[] = { 0,0,1,1, 0,2,1,2, 2 } ;
01161 int n1, n2 ;
01162
01163 int i, j, k, q ;
01164
01165 static double xs[ndm][ndm] ;
01166 static double sx[ndm][ndm] ;
01167
01168 double ss = coor[0] ;
01169 double tt = coor[1] ;
01170 double tempSS, tempTT ;
01171
01172
01173 for ( q=0; q<nNode; q++ ) {
01174
01175 n1 = node1[q] ;
01176 n2 = node2[q] ;
01177
01178 tempSS = this->shape1d(1,n1,ss) ;
01179 tempTT = this->shape1d(1,n2,tt) ;
01180
01181
01182 shp[2][q] = tempSS * tempTT ;
01183
01184
01185 shp[0][q] = (this->shape1d(0,n1,ss)) * tempTT ;
01186
01187
01188 shp[1][q] = tempSS * (this->shape1d(0,n2,tt)) ;
01189
01190 }
01191
01192
01193
01194 for ( i=0; i<ndm; i++ ) {
01195 for ( j=0; j<ndm; j++ ) {
01196
01197 xs[i][j] = 0.0 ;
01198
01199 for ( k=0; k<nNode; k++ )
01200 xs[i][j] += ( x[i][k] * shp[j][k] ) ;
01201
01202 }
01203 }
01204
01205
01206 xsj = xs[0][0]*xs[1][1] - xs[0][1]*xs[1][0] ;
01207
01208 double xsjInv = 1.0/xsj ;
01209
01210 sx[0][0] = xs[1][1] * xsjInv ;
01211 sx[1][1] = xs[0][0] * xsjInv ;
01212 sx[0][1] = -xs[0][1] * xsjInv ;
01213 sx[1][0] = -xs[1][0] * xsjInv ;
01214
01215
01216
01217 double temp ;
01218 for ( i=0; i<nNode; i++ ) {
01219 temp = shp[0][i]*sx[0][0] + shp[1][i]*sx[1][0] ;
01220 shp[1][i] = shp[0][i]*sx[0][1] + shp[1][i]*sx[1][1] ;
01221 shp[0][i] = temp ;
01222 }
01223
01224
01225 return ;
01226 }
01227
01228
01229
01230
01231 double
01232 NineNodeMixedQuad::shape1d( int code, int node, double xi )
01233 {
01234
01235
01236
01237
01238
01239
01240
01241 double result ;
01242
01243 static double oneHalf = 0.50 ;
01244
01245
01246 if ( code == 1 ) {
01247
01248 if ( node == 0 )
01249 result = oneHalf*xi*(xi-1.0) ;
01250
01251 if ( node == 1 )
01252 result = oneHalf*xi*(xi+1.0) ;
01253
01254 if ( node == 2 )
01255 result = 1.0 - xi*xi ;
01256
01257 }
01258
01259 else if ( code == 0 ) {
01260
01261 if ( node == 0 )
01262 result = oneHalf * ( 2.0*xi - 1.0 ) ;
01263
01264 if ( node == 1 )
01265 result = oneHalf * ( 2.0*xi + 1.0 ) ;
01266
01267 if ( node == 2 )
01268 result = -2.0*xi ;
01269
01270 }
01271
01272 return result ;
01273 }
01274
01275
01276
01277 int
01278 NineNodeMixedQuad::displaySelf(Renderer &theViewer, int displayMode, float fact)
01279 {
01280
01281
01282
01283 const Vector &end1Crd = nodePointers[0]->getCrds();
01284 const Vector &end2Crd = nodePointers[4]->getCrds();
01285 const Vector &end3Crd = nodePointers[1]->getCrds();
01286 const Vector &end4Crd = nodePointers[5]->getCrds();
01287 const Vector &end5Crd = nodePointers[2]->getCrds();
01288 const Vector &end6Crd = nodePointers[6]->getCrds();
01289 const Vector &end7Crd = nodePointers[3]->getCrds();
01290 const Vector &end8Crd = nodePointers[7]->getCrds();
01291
01292 const Vector &end1Disp = nodePointers[0]->getDisp();
01293 const Vector &end2Disp = nodePointers[4]->getDisp();
01294 const Vector &end3Disp = nodePointers[1]->getDisp();
01295 const Vector &end4Disp = nodePointers[5]->getDisp();
01296 const Vector &end5Disp = nodePointers[2]->getDisp();
01297 const Vector &end6Disp = nodePointers[6]->getDisp();
01298 const Vector &end7Disp = nodePointers[3]->getDisp();
01299 const Vector &end8Disp = nodePointers[7]->getDisp();
01300
01301
01302 static Matrix coords(8,3) ;
01303 static Vector values(8) ;
01304 static Vector P(8) ;
01305
01306 coords.Zero( ) ;
01307
01308 values(0) = 1 ;
01309 values(1) = 1 ;
01310 values(2) = 1 ;
01311 values(3) = 1 ;
01312 values(4) = 1 ;
01313 values(5) = 1 ;
01314 values(6) = 1 ;
01315 values(7) = 1 ;
01316
01317
01318 if (displayMode < 3 && displayMode > 0)
01319 P = this->getResistingForce();
01320
01321 for (int i = 0; i < 2; i++) {
01322 coords(0,i) = end1Crd(i) + end1Disp(i)*fact;
01323 coords(1,i) = end2Crd(i) + end2Disp(i)*fact;
01324 coords(2,i) = end3Crd(i) + end3Disp(i)*fact;
01325 coords(3,i) = end4Crd(i) + end4Disp(i)*fact;
01326 coords(4,i) = end5Crd(i) + end5Disp(i)*fact;
01327 coords(5,i) = end6Crd(i) + end6Disp(i)*fact;
01328 coords(6,i) = end7Crd(i) + end7Disp(i)*fact;
01329 coords(7,i) = end8Crd(i) + end8Disp(i)*fact;
01330
01331
01332
01333
01334 }
01335
01336
01337 int error = 0;
01338
01339 error += theViewer.drawPolygon (coords, values);
01340
01341 return error;
01342 }
01343
01344
01345
01346
01347 Response*
01348 NineNodeMixedQuad::setResponse(const char **argv, int argc,
01349 Information &eleInfo, OPS_Stream &output)
01350 {
01351 Response *theResponse =0;
01352
01353 output.tag("ElementOutput");
01354 output.attr("eleType","NineNodeMixedQuad");
01355 output.attr("eleTag",this->getTag());
01356 output.attr("node1",connectedExternalNodes[0]);
01357 output.attr("node2",connectedExternalNodes[1]);
01358 output.attr("node3",connectedExternalNodes[2]);
01359 output.attr("node4",connectedExternalNodes[3]);
01360 output.attr("node5",connectedExternalNodes[4]);
01361 output.attr("node6",connectedExternalNodes[5]);
01362 output.attr("node7",connectedExternalNodes[6]);
01363 output.attr("node8",connectedExternalNodes[7]);
01364 output.attr("node9",connectedExternalNodes[8]);
01365
01366 char dataOut[10];
01367 if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0) {
01368
01369 for (int i=1; i<=4; i++) {
01370 sprintf(dataOut,"P1_%d",i);
01371 output.tag("ResponseType",dataOut);
01372 sprintf(dataOut,"P2_%d",i);
01373 output.tag("ResponseType",dataOut);
01374 }
01375
01376 theResponse = new ElementResponse(this, 1, resid);
01377 } else if (strcmp(argv[0],"material") == 0 || strcmp(argv[0],"integrPoint") == 0) {
01378 int pointNum = atoi(argv[1]);
01379 if (pointNum > 0 && pointNum <= 9) {
01380
01381 output.tag("GaussPoint");
01382 output.attr("number",pointNum);
01383 output.attr("eta",sg[pointNum-1]);
01384 output.attr("neta",sg[pointNum-1]);
01385
01386 theResponse = materialPointers[pointNum-1]->setResponse(&argv[2], argc-2, eleInfo, output);
01387
01388 output.endTag();
01389
01390 } else if (strcmp(argv[0],"stresses") ==0) {
01391
01392 for (int i=0; i<9; i++) {
01393 output.tag("GaussPoint");
01394 output.attr("number",i+1);
01395 output.attr("eta",sg[i]);
01396 output.attr("neta",sg[i]);
01397
01398 output.tag("NdMaterialOutput");
01399 output.attr("classType", materialPointers[i]->getClassTag());
01400 output.attr("tag", materialPointers[i]->getTag());
01401
01402 output.tag("ResponseType","UnknownStress");
01403 output.tag("ResponseType","UnknownStress");
01404 output.tag("ResponseType","UnknownStress");
01405 output.tag("ResponseType","UnknownStress");
01406
01407 output.endTag();
01408 output.endTag();
01409 }
01410
01411 theResponse = new ElementResponse(this, 3, Vector(4*9));
01412 }
01413 }
01414
01415 output.endTag();
01416
01417 return theResponse;
01418 }
01419
01420 int
01421 NineNodeMixedQuad::getResponse(int responseID, Information &eleInfo)
01422 {
01423 if (responseID == 1) {
01424
01425 return eleInfo.setVector(this->getResistingForce());
01426
01427 } else if (responseID == 3) {
01428
01429
01430 static Vector stresses(4*9);
01431 int cnt = 0;
01432 for (int i = 0; i < 9; i++) {
01433
01434
01435 const Vector &sigma = materialPointers[i]->getStress();
01436 stresses(cnt) = sigma(0);
01437 stresses(cnt+1) = sigma(1);
01438 stresses(cnt+2) = sigma(2);
01439 stresses(cnt+3) = sigma(3);
01440 cnt += 4;
01441 }
01442 return eleInfo.setVector(stresses);
01443
01444 } else
01445
01446 return -1;
01447 }
01448
01449
01450 int
01451 NineNodeMixedQuad::sendSelf (int commitTag, Channel &theChannel)
01452 {
01453 int res = 0;
01454
01455
01456
01457
01458 int dataTag = this->getDbTag();
01459
01460
01461
01462 int matDbTag;
01463
01464 static ID idData(28);
01465
01466 int i;
01467 for (i = 0; i < 9; i++) {
01468 idData(i) = materialPointers[i]->getClassTag();
01469 matDbTag = materialPointers[i]->getDbTag();
01470
01471
01472 if (matDbTag == 0) {
01473 matDbTag = theChannel.getDbTag();
01474 if (matDbTag != 0)
01475 materialPointers[i]->setDbTag(matDbTag);
01476 }
01477 idData(i+9) = matDbTag;
01478 }
01479
01480 idData(18) = this->getTag();
01481 idData(19) = connectedExternalNodes(0);
01482 idData(20) = connectedExternalNodes(1);
01483 idData(21) = connectedExternalNodes(2);
01484 idData(22) = connectedExternalNodes(3);
01485 idData(23) = connectedExternalNodes(4);
01486 idData(24) = connectedExternalNodes(5);
01487 idData(25) = connectedExternalNodes(6);
01488 idData(26) = connectedExternalNodes(7);
01489 idData(27) = connectedExternalNodes(8);
01490
01491
01492 res += theChannel.sendID(dataTag, commitTag, idData);
01493 if (res < 0) {
01494 opserr << "WARNING NineNodeMixedQuad::sendSelf() - " << this->getTag() << " failed to send ID\n";
01495
01496 return res;
01497 }
01498
01499
01500 for (i = 0; i < 9; i++) {
01501 res += materialPointers[i]->sendSelf(commitTag, theChannel);
01502 if (res < 0) {
01503 opserr << "WARNING NineNodeMixedQuad::sendSelf() - " << this->getTag() << " failed to send its Material\n";
01504 return res;
01505 }
01506 }
01507
01508 return res;
01509 }
01510
01511 int
01512 NineNodeMixedQuad::recvSelf (int commitTag,
01513 Channel &theChannel,
01514 FEM_ObjectBroker &theBroker)
01515 {
01516 int res = 0;
01517
01518 int dataTag = this->getDbTag();
01519
01520 static ID idData(28);
01521
01522 res += theChannel.recvID(dataTag, commitTag, idData);
01523 if (res < 0) {
01524 opserr << "WARNING NineNodeMixedQuad::recvSelf() - " << this->getTag() << " failed to receive ID\n";
01525 return res;
01526 }
01527
01528 this->setTag(idData(18));
01529 connectedExternalNodes(0) = idData(19);
01530 connectedExternalNodes(1) = idData(20);
01531 connectedExternalNodes(2) = idData(21);
01532 connectedExternalNodes(3) = idData(22);
01533 connectedExternalNodes(4) = idData(23);
01534 connectedExternalNodes(5) = idData(24);
01535 connectedExternalNodes(6) = idData(25);
01536 connectedExternalNodes(7) = idData(26);
01537 connectedExternalNodes(8) = idData(27);
01538
01539
01540 int i;
01541
01542 if (materialPointers[0] == 0) {
01543 for (i = 0; i < 9; i++) {
01544 int matClassTag = idData(i);
01545 int matDbTag = idData(i+9);
01546
01547 materialPointers[i] = theBroker.getNewNDMaterial(matClassTag);
01548 if (materialPointers[i] == 0) {
01549 opserr << "NineNodeMixedQuad::recvSelf() - Broker could not create NDMaterial of class type" << matClassTag << endln;
01550 return -1;
01551 }
01552
01553 materialPointers[i]->setDbTag(matDbTag);
01554 res += materialPointers[i]->recvSelf(commitTag, theChannel, theBroker);
01555 if (res < 0) {
01556 opserr << "NineNodeMixedQuad::recvSelf() - material " <<
01557 i << "failed to recv itself\n";
01558 return res;
01559 }
01560 }
01561 }
01562
01563 else {
01564 for (i = 0; i < 9; i++) {
01565 int matClassTag = idData(i);
01566 int matDbTag = idData(i+9);
01567
01568
01569 if (materialPointers[i]->getClassTag() != matClassTag) {
01570 delete materialPointers[i];
01571 materialPointers[i] = theBroker.getNewNDMaterial(matClassTag);
01572 if (materialPointers[i] == 0) {
01573 opserr << "NineNodeMixedQuad::recvSelf() - Broker could not create NDMaterial of class type" << matClassTag << endln;
01574 exit(-1);
01575 }
01576 }
01577
01578 materialPointers[i]->setDbTag(matDbTag);
01579 res += materialPointers[i]->recvSelf(commitTag, theChannel, theBroker);
01580 if (res < 0) {
01581 opserr << "NineNodeMixedQuad::recvSelf() - material " <<
01582 i << "failed to recv itself\n";
01583 return res;
01584 }
01585 }
01586 }
01587
01588 return res;
01589 }
01590