NineNodeMixedQuad.cpp

Go to the documentation of this file.
00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.11 $
00022 // $Date: 2006/08/04 19:07:15 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/fourNodeQuad/NineNodeMixedQuad.cpp,v $
00024 
00025 // Ed "C++" Love
00026 //
00027 // Mixed Presssure/Volume Nine Node Quadrilateral
00028 // Plane Strain (NOT PLANE STRESS)
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 //static data
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 //quadrature data
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 //null constructor
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 //full constructor
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     } //end if
00107       
00108   } //end for i 
00109 }
00110 
00111 
00112 //destructor 
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   } //end for i
00123 
00124   if (load != 0)
00125     delete load;
00126 
00127   if (Ki != 0)
00128     delete Ki;
00129 }
00130 
00131 
00132 //set domain
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 //get the number of external nodes
00144 int 
00145 NineNodeMixedQuad::getNumExternalNodes( ) const
00146 {
00147   return 9 ;
00148 } 
00149  
00150 
00151 //return connected external nodes
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 //return number of dofs
00166 int 
00167 NineNodeMixedQuad::getNumDOF( ) 
00168 {
00169   return 18 ;
00170 }
00171 
00172 
00173 //commit state
00174 int 
00175 NineNodeMixedQuad::commitState( )
00176 {
00177   int i ;
00178   int success = 0 ;
00179 
00180   // call element commitState to do any base class stuff
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 //revert to last commit 
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 //revert to start 
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 //print out element data
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 //return stiffness matrix 
00246 const Matrix& 
00247 NineNodeMixedQuad::getTangentStiff( ) 
00248 {
00249   int tang_flag = 1 ; //get the tangent 
00250 
00251   //do tangent and residual here
00252   formResidAndTangent( tang_flag ) ;  
00253 
00254   return stiff ;
00255 }    
00256 
00257 
00258 //return secant matrix 
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 ;  // determinant jacaobian matrix 
00285 
00286   static double dvol[numberGauss] ; //volume element
00287 
00288   static double gaussPoint[ndm] ;
00289 
00290   static double natCoorArray[ndm][numberGauss] ;
00291 
00292   static Vector strain(nstress) ;  //strain
00293 
00294   static double shp[nShape][numberNodes] ;  //shape functions at a gauss point
00295 
00296   static double Shape[nShape][numberNodes][numberGauss] ; //all the shape functions
00297 
00298   static double shpBar[nShape][numberNodes][nMixed] ; //mean value of shape functions
00299 
00300   static double rightHandSide[nShape][numberNodes][nMixed] ;
00301 
00302   static Vector residJ(ndf) ; //nodeJ residual 
00303 
00304   static Matrix stiffJK(ndf,ndf) ; //nodeJK stiffness 
00305 
00306   static Vector stress(nstress) ;  //stress
00307 
00308   static Matrix dd(nstress,nstress) ;  //material tangent
00309 
00310   static double interp[nMixed] ;
00311 
00312   static Matrix Proj(3,3) ;   //projection matrix 
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   //---------B-matrices------------------------------------
00321 
00322     static Matrix BJ(nstress,ndf) ;      // B matrix node J
00323 
00324     static Matrix BJtran(ndf,nstress) ;
00325 
00326     static Matrix BK(nstress,ndf) ;      // B matrix node k
00327 
00328     static Matrix BJtranD(ndf,nstress) ;
00329 
00330   //-------------------------------------------------------
00331 
00332   
00333   //zero stiffness and residual 
00334   stiff.Zero( ) ;
00335 
00336   //node coordinates
00337   computeBasis() ;
00338 
00339   //zero mean shape functions
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     }//end for q
00349   } // end for p
00350 
00351 
00352   //zero volume
00353   volume = 0.0 ;
00354 
00355   //zero projection matrix  
00356   Proj.Zero( ) ;
00357   ProjInv.Zero( ) ;
00358 
00359   //gauss loop to compute and save shape functions 
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         //save gauss point locations
00370         natCoorArray[0][count] = gaussPoint[0] ;
00371         natCoorArray[1][count] = gaussPoint[1] ;
00372 
00373 
00374         //get shape functions    
00375         shape2dNine( gaussPoint, xl, shp, xsj ) ;
00376 
00377 
00378         //save shape functions
00379         for ( p=0; p<nShape; p++ ) {
00380           for ( q=0; q<numberNodes; q++ )
00381             Shape[p][q][count] = shp[p][q] ;
00382         } // end for p
00383 
00384         
00385         //volume element to also be saved
00386         dvol[count] = ( wg[i]*wg[j] ) * xsj ;  
00387 
00388 
00389         //add to projection matrix
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         }//end for r
00398 
00399         volume += dvol[count] ;
00400         
00401         
00402         //add to mean shape functions
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           }//end for q 
00410         } // end for p
00411 
00412 
00413         //increment gauss point counter
00414         count++ ;
00415 
00416     } //end for j
00417   } // end for i 
00418   
00419 
00420 
00421   //invert projection matrix
00422   //int Solve(const Matrix &M, Matrix &res) const;
00423   Proj.Solve( Iden, ProjInv ) ;
00424   
00425   //mean value of shape functions
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       }//end for r
00433 
00434     }//end for q
00435   }//end for p
00436 
00437 
00438   //gauss loop 
00439   for ( i=0; i<numberGauss; i++ ) {
00440     
00441     //extract gauss point location
00442     gaussPoint[0] = natCoorArray[0][i] ;
00443     gaussPoint[1] = natCoorArray[1][i] ;
00444 
00445     //extract shape functions from saved array
00446     for ( p=0; p<nShape; p++ ) {
00447        for ( q=0; q<numberNodes; q++ )
00448           shp[p][q]  = Shape[p][q][i] ;
00449     } // end for p
00450 
00451     dd = materialPointers[i]->getInitialTangent( ) ;
00452     dd *= dvol[i] ;
00453 
00454 
00455     //residual and tangent calculations node loops
00456 
00457     jj = 0 ;
00458     for ( j=0; j<numberNodes; j++ ) {
00459 
00460       BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00461    
00462       //transpose 
00463       //BJtran = transpose( nstress, ndf, BJ ) ;
00464       for (p=0; p<ndf; p++) {
00465         for (q=0; q<nstress; q++) 
00466           BJtran(p,q) = BJ(q,p) ;
00467       }//end for p
00468 
00469 
00470       //BJtranD = BJtran * dd ;
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         //stiffJK =  BJtranD * BK  ;
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         } //end for p
00486         
00487         kk += ndf ;
00488       }//end for k loop
00489 
00490       jj += ndf ;
00491     }//end for j loop
00492   }//end for i gauss loop 
00493 
00494   Ki = new Matrix(stiff);
00495 
00496   return stiff;
00497 }
00498     
00499 
00500 //return mass matrix
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   // check to see if have mass
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   // Compute mass matrix
00549   int tangFlag = 1 ;
00550   formInertiaTerms( tangFlag ) ;
00551 
00552   // store computed RV fro nodes in resid vector
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   // create the load vector if one does not exist
00562   if (load == 0) 
00563     load = new Vector(numberNodes*ndf);
00564 
00565   // add -M * RV(accel) to the load vector
00566   load->addMatrixVector(1.0, mass, resid, -1.0);
00567   
00568   return 0;
00569 }
00570 
00571 
00572 //get residual
00573 const Vector& 
00574 NineNodeMixedQuad::getResistingForce( ) 
00575 {
00576   int tang_flag = 0 ; //don't get the tangent
00577 
00578   formResidAndTangent( tang_flag ) ;
00579 
00580   // subtract external loads 
00581   if (load != 0)
00582     resid -= *load;
00583 
00584   return resid ;   
00585 }
00586 
00587 
00588 //get residual with inertia terms
00589 const Vector& 
00590 NineNodeMixedQuad::getResistingForceIncInertia( )
00591 {
00592   int tang_flag = 0 ; //don't get the tangent
00593 
00594   static Vector res(18);
00595 
00596   //do tangent and residual here 
00597   formResidAndTangent( tang_flag ) ;
00598 
00599   //inertia terms
00600   formInertiaTerms( tang_flag ) ;
00601 
00602   res = resid;
00603 
00604   // add the damping forces if rayleigh damping
00605   if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
00606     res += this->getRayleighDampingForces();
00607 
00608   // subtract external loads 
00609   if (load != 0)
00610     res -= *load;
00611 
00612   return res;
00613 }
00614 
00615 //*****************************************************************************
00616 //form inertia terms
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   //  static const int numberGauss = 9 ;
00629 
00630   static const int nShape = 3 ;
00631 
00632   static const int massIndex = nShape - 1 ;
00633 
00634   double xsj ;  // determinant jacaobian matrix 
00635 
00636   double dvol ; //volume element
00637 
00638   static double shp[nShape][numberNodes] ;  //shape functions at a gauss point
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   //zero mass 
00653   mass.Zero( ) ;
00654 
00655   //node coordinates
00656   computeBasis() ;
00657 
00658   //gauss loop 
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     //get shape functions    
00667     shape2dNine( GaussPoint, xl, shp, xsj ) ;
00668 
00669     //volume element
00670     dvol = ( wg[p] * wg[q] ) * xsj ;
00671 
00672     //node loop to compute acceleration
00673     momentum.Zero( ) ;
00674     for ( j = 0; j < numberNodes; j++ ) 
00675       //momentum += shp[massIndex][j] * ( nodePointers[j]->getTrialAccel()  ) ; 
00676       momentum.addVector( 1.0,
00677                           nodePointers[j]->getTrialAccel(),
00678                           shp[massIndex][j] ) ;
00679 
00680 
00681     //density
00682     rho = materialPointers[count]->getRho() ;
00683 
00684     //multiply acceleration by density to form momentum
00685     momentum *= rho ;
00686 
00687 
00688     //residual and tangent calculations node loops
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          //multiply by density
00700          temp *= rho ;
00701 
00702          //node-node mass
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           } // end for k loop
00711 
00712       } // end if tang_flag 
00713 
00714     } // end for j loop
00715 
00716     
00717     count++ ;
00718     }//end for q gauss loop
00719   } //end for p gauss loop 
00720 
00721 }
00722 
00723 //*********************************************************************
00724 //form residual and tangent
00725 
00726 void  
00727 NineNodeMixedQuad::formResidAndTangent( int tang_flag ) 
00728 {
00729 
00730   //strains ordered : eps11, eps22, eps33, 2*eps12 
00731   //volumtric strains projected onto {1, \xi, \eta} natural coordinates
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 ;  // determinant jacaobian matrix 
00755 
00756   static double dvol[numberGauss] ; //volume element
00757 
00758   static double gaussPoint[ndm] ;
00759 
00760   static double natCoorArray[ndm][numberGauss] ;
00761 
00762   static Vector strain(nstress) ;  //strain
00763 
00764   static double shp[nShape][numberNodes] ;  //shape functions at a gauss point
00765 
00766   static double Shape[nShape][numberNodes][numberGauss] ; //all the shape functions
00767 
00768   static double shpBar[nShape][numberNodes][nMixed] ; //mean value of shape functions
00769 
00770   static double rightHandSide[nShape][numberNodes][nMixed] ;
00771 
00772   static Vector residJ(ndf) ; //nodeJ residual 
00773 
00774   static Matrix stiffJK(ndf,ndf) ; //nodeJK stiffness 
00775 
00776   static Vector stress(nstress) ;  //stress
00777 
00778   static Matrix dd(nstress,nstress) ;  //material tangent
00779 
00780   static double interp[nMixed] ;
00781 
00782   static Matrix Proj(3,3) ;   //projection matrix 
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   //---------B-matrices------------------------------------
00791 
00792     static Matrix BJ(nstress,ndf) ;      // B matrix node J
00793 
00794     static Matrix BJtran(ndf,nstress) ;
00795 
00796     static Matrix BK(nstress,ndf) ;      // B matrix node k
00797 
00798     static Matrix BJtranD(ndf,nstress) ;
00799 
00800   //-------------------------------------------------------
00801 
00802   
00803   //zero stiffness and residual 
00804   stiff.Zero( ) ;
00805   resid.Zero( ) ;
00806 
00807   //node coordinates
00808   computeBasis() ;
00809 
00810   //zero mean shape functions
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     }//end for q
00820   } // end for p
00821 
00822 
00823   //zero volume
00824   volume = 0.0 ;
00825 
00826   //zero projection matrix  
00827   Proj.Zero( ) ;
00828   ProjInv.Zero( ) ;
00829 
00830 
00831   //gauss loop to compute and save shape functions 
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         //save gauss point locations
00842         natCoorArray[0][count] = gaussPoint[0] ;
00843         natCoorArray[1][count] = gaussPoint[1] ;
00844 
00845 
00846         //get shape functions    
00847         shape2dNine( gaussPoint, xl, shp, xsj ) ;
00848 
00849 
00850         //save shape functions
00851         for ( p=0; p<nShape; p++ ) {
00852           for ( q=0; q<numberNodes; q++ )
00853             Shape[p][q][count] = shp[p][q] ;
00854         } // end for p
00855 
00856         
00857         //volume element to also be saved
00858         dvol[count] = ( wg[i]*wg[j] ) * xsj ;  
00859 
00860 
00861         //add to projection matrix
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         }//end for r
00870 
00871         volume += dvol[count] ;
00872         
00873         
00874         //add to mean shape functions
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           }//end for q 
00882         } // end for p
00883 
00884 
00885         //increment gauss point counter
00886         count++ ;
00887 
00888     } //end for j
00889   } // end for i 
00890   
00891 
00892 
00893   //invert projection matrix
00894   //int Solve(const Matrix &M, Matrix &res) const;
00895   Proj.Solve( Iden, ProjInv ) ;
00896   
00897   //mean value of shape functions
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       }//end for r
00905 
00906     }//end for q
00907   }//end for p
00908 
00909 
00910   //gauss loop 
00911   for ( i=0; i<numberGauss; i++ ) {
00912     
00913     //extract gauss point location
00914     gaussPoint[0] = natCoorArray[0][i] ;
00915     gaussPoint[1] = natCoorArray[1][i] ;
00916 
00917     //extract shape functions from saved array
00918     for ( p=0; p<nShape; p++ ) {
00919        for ( q=0; q<numberNodes; q++ )
00920           shp[p][q]  = Shape[p][q][i] ;
00921     } // end for p
00922 
00923 
00924     //zero the strains
00925     strain.Zero( ) ;
00926 
00927     // j-node loop to compute strain 
00928     for ( j=0; j<numberNodes; j++ )  {
00929 
00930       //compute B matrix 
00931 
00932       BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00933       
00934       //nodal displacements 
00935       const Vector &ul = nodePointers[j]->getTrialDisp( ) ;
00936 
00937       //compute the strain
00938       //strain += (BJ*ul) ; 
00939       strain.addMatrixVector(1.0,  BJ,ul,1.0 ) ;
00940 
00941     } // end for j
00942   
00943 
00944 
00945     //send the strain to the material 
00946     success = materialPointers[i]->setTrialStrain( strain ) ;
00947 
00948     //compute the stress
00949     stress = materialPointers[i]->getStress( ) ;
00950 
00951 
00952     //multiply by volume element
00953     stress  *= dvol[i] ;
00954 
00955     if ( tang_flag == 1 ) {
00956       dd = materialPointers[i]->getTangent( ) ;
00957       dd *= dvol[i] ;
00958     } //end if tang_flag
00959 
00960 
00961     //residual and tangent calculations node loops
00962 
00963     jj = 0 ;
00964     for ( j=0; j<numberNodes; j++ ) {
00965 
00966       BJ = computeBbar( j, gaussPoint, shp, shpBar ) ;
00967    
00968       //transpose 
00969       //BJtran = transpose( nstress, ndf, BJ ) ;
00970       for (p=0; p<ndf; p++) {
00971         for (q=0; q<nstress; q++) 
00972           BJtran(p,q) = BJ(q,p) ;
00973       }//end for p
00974 
00975 
00976       //residual
00977       //residJ = BJtran * stress ;
00978       residJ.addMatrixVector(0.0,  BJtran,stress,1.0);
00979 
00980       //residual 
00981       for ( p=0; p<ndf; p++ )
00982         resid( jj + p ) += residJ(p)  ;
00983 
00984 
00985       if ( tang_flag == 1 ) {
00986 
00987         //BJtranD = BJtran * dd ;
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             //stiffJK =  BJtranD * BK  ;
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             } //end for p
01003 
01004             kk += ndf ;
01005          }//end for k loop
01006 
01007       }//end if tang_flag 
01008 
01009       jj += ndf ;
01010     }//end for j loop
01011 
01012 
01013   }//end for i gauss loop 
01014 
01015   
01016   return ;
01017 }
01018 
01019 
01020 //*************************************************************************
01021 //compute local coordinates and basis
01022 
01023 void   
01024 NineNodeMixedQuad::computeBasis( ) 
01025 {
01026 
01027   //nodal coordinates 
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   }  //end for i 
01037 
01038 }
01039 
01040 //*************************************************************************
01041 //compute Bbar
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 //---B Matrices in standard {1,2,3} mechanics notation---------
01065 //
01066 //                -                  -
01067 //               |  2N,1    -N,2      | 
01068 // Bdev =  (1/3) |  -N,1    2N,2      |  (3x2)
01069 //               |  -N,1    -N,2      |   
01070 //                -                  -
01071 //
01072 //                -                 -
01073 //               |  N,1      N,2     | 
01074 // Bvol =  (1/3) |  N,1      N,2     |  (3x2)
01075 //               |  N,1      N,2     |   
01076 //                -                 -
01077 //
01078 //                -                   -
01079 //               |                     |
01080 //               |    Bdev + Bvol      |
01081 //   B       =   |                     | 
01082 //               |---------------------|   (4x2)
01083 //               |   N,2     N,1       |
01084 //                -                   -       
01085 //
01086 //---------------------------------------------------------------
01087 
01088 
01089   Bbar.Zero( ) ;
01090 
01091   //deviatoric
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   //volume interpolation functions
01103   interp[0] = 1.0 ;
01104   interp[1] = natCoor[0] ;
01105   interp[2] = natCoor[1] ; 
01106 
01107   //volumetric 
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   }//end for i
01117 
01118   //standard displacement formulation
01119   //c0 = shp[0][node] ;
01120   //c1 = shp[1][node] ;
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   //extensional terms
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   }//end for i
01137 
01138 
01139   //shear terms
01140   Bbar(3,0) = shp[1][node] ;
01141   Bbar(3,1) = shp[0][node] ;
01142 
01143 
01144   return Bbar ;
01145 }
01146 
01147 //**************************************************************************
01148 //shape function routine for nine node quads
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   //shape functions and derivatives in natural coordinates
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     //function
01182     shp[2][q] =  tempSS * tempTT ;
01183 
01184     //ss-derivative
01185     shp[0][q] =  (this->shape1d(0,n1,ss)) * tempTT ;
01186     
01187     //tt-derivative
01188     shp[1][q] =  tempSS * (this->shape1d(0,n2,tt)) ;  
01189 
01190   }//end for q
01191 
01192 
01193   // Construct jacobian and its inverse
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     } //end for j
01203   }  // end for i 
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   //form global derivatives 
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   } // end for i
01223 
01224 
01225   return ;
01226 }
01227         
01228 //***********************************************************************
01229 //1d quadratic shape functions
01230 
01231 double 
01232 NineNodeMixedQuad::shape1d( int code, int node, double xi ) 
01233 {
01234 
01235   //one-dimensional quadratic shape functions
01236   //
01237   // o------o------o
01238   // 0      2      1
01239   //
01240 
01241   double result ;
01242 
01243   static double oneHalf = 0.50 ;
01244 
01245   //shape functions
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   //shape function derivatives
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   }//end if
01271 
01272   return result ;
01273 }
01274 
01275 //***********************************************************************
01276 
01277 int
01278 NineNodeMixedQuad::displaySelf(Renderer &theViewer, int displayMode, float fact)
01279 {
01280     // first determine the end points of the quad based on
01281     // the display factor (a measure of the distorted image)
01282     // store this information in 4 3d vectors v1 through v4
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       /*      if (displayMode < 3 && displayMode > 0)
01331         values(i) = P(displayMode*2+i);
01332       else
01333       values(i) = 1;  */
01334     }
01335 
01336     //opserr << coords;
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(); // GaussPoint
01408         output.endTag(); // NdMaterialOutput
01409       }
01410 
01411       theResponse =  new ElementResponse(this, 3, Vector(4*9));
01412     }
01413   }
01414         
01415   output.endTag(); // ElementOutput
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     // Loop over the integration points
01430     static Vector stresses(4*9);
01431     int cnt = 0;
01432     for (int i = 0; i < 9; i++) {
01433 
01434       // Get material stress response
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   // note: we don't check for dataTag == 0 for Element
01456   // objects as that is taken care of in a commit by the Domain
01457   // object - don't want to have to do the check if sending data
01458   int dataTag = this->getDbTag();
01459   
01460 
01461   // Now quad sends the ids of its materials
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     // NOTE: we do have to ensure that the material has a database
01471     // tag if we are sending to a database channel.
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   // Finally, quad asks its material objects to send themselves
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   // Quad now receives the tags of its four external nodes
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       // Allocate new material with the sent class tag
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       // Now receive materials into the newly allocated space
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   // Number of materials is the same, receive materials into current space
01563   else {
01564     for (i = 0; i < 9; i++) {
01565       int matClassTag = idData(i);
01566       int matDbTag = idData(i+9);
01567       // Check that material is of the right type; if not,
01568       // delete it and create a new one of the right type
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       // Receive the material
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 

Generated on Mon Oct 23 15:05:07 2006 for OpenSees by doxygen 1.5.0