OpenSees » SRC » element » UP-ucsd

Nine_Four_Node_QuadUPOld.cpp

Go to the documentation of this file.
00001 
00002 // Description: This file contains the class definition for                  //
00003 // NineFourNodeQuadUP, a 9-4-node (9 node for solid and 4 node for fluid) //
00004 // plane strain element for solid-fluid fully coupled analysis. This         //
00005 // implementation is a simplified u-p formulation of Biot theory             //
00006 // (u - solid displacement, p - fluid pressure). Each element node has two   //
00007 // DOFs for u and 1 DOF for p.                                               //
00008 //                                                                           //
00009 // Written by Zhaohui Yang      (March 2004)                                     //
00010 //                                                                           //
00012 
00013 // $Revision: 1.1 $
00014 // $Date: 2005/09/22 21:28:36 $
00015 // $Source: /usr/local/cvs/OpenSees/SRC/element/UP-ucsd/Nine_Four_Node_QuadUPOld.cpp,v $
00016 
00017 #include <Nine_Four_Node_QuadUP.h>
00018 #include <Node.h>
00019 #include <NDMaterial.h>
00020 #include <Matrix.h>
00021 #include <Vector.h>
00022 #include <ID.h>
00023 #include <Renderer.h>
00024 #include <Domain.h>
00025 #include <string.h>
00026 #include <Information.h>
00027 #include <Channel.h>
00028 #include <FEM_ObjectBroker.h>
00029 #include <ElementResponse.h>
00030 
00031 Matrix NineFourNodeQuadUP::K(22,22);
00032 Vector NineFourNodeQuadUP::P(22);
00033 double NineFourNodeQuadUP::shgu[3][9][9];       
00034 double NineFourNodeQuadUP::shgp[3][4][4];       
00035 double NineFourNodeQuadUP::shgq[3][9][4];       
00036 double NineFourNodeQuadUP::shlu[3][9][9];       
00037 double NineFourNodeQuadUP::shlp[3][4][4];       
00038 double NineFourNodeQuadUP::shlq[3][9][4];       
00039 double NineFourNodeQuadUP::wu[9];       
00040 double NineFourNodeQuadUP::wp[4];       
00041 double NineFourNodeQuadUP::dvolu[9];  
00042 double NineFourNodeQuadUP::dvolp[4];  
00043 double NineFourNodeQuadUP::dvolq[4];  
00044 const int NineFourNodeQuadUP::nintu=9;
00045 const int NineFourNodeQuadUP::nintp=4;
00046 const int NineFourNodeQuadUP::nenu=9;
00047 const int NineFourNodeQuadUP::nenp=4;
00048  
00049 NineFourNodeQuadUP::NineFourNodeQuadUP(int tag, 
00050         int nd1, int nd2, int nd3, int nd4,int nd5, int nd6, int nd7, int nd8,int nd9,
00051         NDMaterial &m, const char *type, double t, double bulk, double r,
00052                   double p1, double p2, double b1, double b2)
00053 :Element (tag, ELE_TAG_Nine_Four_Node_QuadUP), 
00054   theMaterial(0), connectedExternalNodes(9), 
00055   Ki(0), Q(22), thickness(t), kc(bulk), rho(r)
00056 {
00057     this->shapeFunction(wu, nintu, nenu, 0);
00058 /*      for( int L = 0; L < nintu; L++) {
00059                 for( int j = 0; j < nenu; j++) {
00060                 printf("%5d %5d %15.6e %15.6e %15.6e\n", L+1, j+1,
00061                         shlu[0][j][L],shlu[1][j][L],shlu[2][j][L]);
00062                 }
00063         }
00064         exit(-1);
00065 */
00066     this->shapeFunction(wp, nintp, nenp, 1);
00067     this->shapeFunction(wp, nintp, nenu, 2);
00068 
00069         // Body forces
00070         b[0] = b1;
00071         b[1] = b2;
00072         // Permeabilities
00073     perm[0] = p1;
00074     perm[1] = p2;
00075 
00076     // Allocate arrays of pointers to NDMaterials
00077     theMaterial = new NDMaterial *[nintu];
00078     
00079     if (theMaterial == 0) {
00080       opserr << "NineFourNodeQuadUP::NineFourNodeQuadUP - failed allocate material model pointer\n";
00081       exit(-1);
00082     }
00083 
00084     for (int i = 0; i < nintu; i++) {
00085       
00086       // Get copies of the material model for each integration point
00087       theMaterial[i] = m.getCopy(type);
00088       
00089       // Check allocation
00090       if (theMaterial[i] == 0) {
00091              opserr << "NineFourNodeQuadUP::NineFourNodeQuadUP -- failed to get a copy of material model\n";
00092              exit(-1);
00093       }
00094     }
00095 
00096     // Set connected external node IDs
00097     connectedExternalNodes(0) = nd1;
00098     connectedExternalNodes(1) = nd2;
00099     connectedExternalNodes(2) = nd3;
00100     connectedExternalNodes(3) = nd4;
00101     connectedExternalNodes(4) = nd5;
00102     connectedExternalNodes(5) = nd6;
00103     connectedExternalNodes(6) = nd7;
00104     connectedExternalNodes(7) = nd8;
00105     connectedExternalNodes(8) = nd9;
00106 }
00107  
00108 
00109 NineFourNodeQuadUP::NineFourNodeQuadUP()
00110 :Element (0,ELE_TAG_Nine_Four_Node_QuadUP),  
00111   theMaterial(0), connectedExternalNodes(9), 
00112   Ki(0), Q(22), thickness(0.0), kc(0.0), rho(0.0)
00113 {
00114     this->shapeFunction(wu, nintu, nenu, 0);
00115     this->shapeFunction(wp, nintp, nenp, 1);
00116     this->shapeFunction(wp, nintp, nenu, 2);
00117 }
00118 
00119 NineFourNodeQuadUP::~NineFourNodeQuadUP()
00120 {    
00121     for (int i = 0; i < nintu; i++) {
00122       if (theMaterial[i])
00123         delete theMaterial[i];
00124     }
00125 
00126     // Delete the array of pointers to NDMaterial pointer arrays
00127     if (theMaterial)
00128                 delete [] theMaterial;
00129 
00130     for (int ii = 0; ii < nenu; ii++) theNodes[ii] = 0 ;
00131 
00132     if (Ki != 0)
00133       delete Ki;
00134 }
00135 
00136 int
00137 NineFourNodeQuadUP::getNumExternalNodes() const
00138 {
00139     return nenu;
00140 }
00141 
00142 const ID&
00143 NineFourNodeQuadUP::getExternalNodes()
00144 {
00145     return connectedExternalNodes;
00146 }
00147 
00148 Node **
00149 NineFourNodeQuadUP::getNodePtrs()
00150 {
00151   return theNodes;
00152 }
00153 
00154 int
00155 NineFourNodeQuadUP::getNumDOF()
00156 {
00157     return 22;
00158 }
00159 
00160 void
00161 NineFourNodeQuadUP::setDomain(Domain *theDomain)
00162 {
00163   // Check Domain is not null - invoked when object removed from a domain
00164   if (theDomain == 0) {
00165     for (int i=0; i<nenu; i++) theNodes[i] = 0;
00166     return;
00167   }
00168   
00169   int i;
00170   for (i=0; i<nenu; i++) {
00171     theNodes[i] = theDomain->getNode(connectedExternalNodes(i));
00172     if (theNodes[i] == 0) {
00173       opserr << "FATAL ERROR NineFourNodeQuadUP, node not found in domain, tag "
00174              << this->getTag();
00175       return;
00176     }
00177   }
00178   
00179   int dof;
00180   for (i=0; i<nenu; i++) {
00181     dof = theNodes[i]->getNumberDOF();
00182     if ((i<nenp && dof != 3) || (i>=nenp && dof != 2)) {
00183       opserr << "FATAL ERROR NineFourNodeQuadUP, has wrong number of DOFs at its nodes "
00184              << this->getTag();
00185       return;
00186                  }
00187   }
00188   
00189   this->DomainComponent::setDomain(theDomain);
00190 }
00191 
00192 
00193 int
00194 NineFourNodeQuadUP::commitState()
00195 {
00196     int retVal = 0;
00197 
00198     // call element commitState to do any base class stuff
00199     if ((retVal = this->Element::commitState()) != 0) {
00200       opserr << "Nine_Four_Node_Quad_UP::commitState () - failed in base class";
00201     }    
00202 
00203     // Loop over the integration points and commit the material states
00204     for (int i = 0; i < nintu; i++)
00205       retVal += theMaterial[i]->commitState();
00206 
00207     return retVal;
00208 }
00209 
00210 int
00211 NineFourNodeQuadUP::revertToLastCommit()
00212 {
00213     int retVal = 0;
00214 
00215     // Loop over the integration points and revert to last committed state
00216     for (int i = 0; i < nintu; i++)
00217                 retVal += theMaterial[i]->revertToLastCommit();
00218 
00219     return retVal;
00220 }
00221 
00222 int
00223 NineFourNodeQuadUP::revertToStart()
00224 {
00225     int retVal = 0;
00226 
00227     // Loop over the integration points and revert states to start
00228     for (int i = 0; i < nintu; i++)
00229                 retVal += theMaterial[i]->revertToStart();
00230 
00231     return retVal;
00232 }
00233 
00234 int
00235 NineFourNodeQuadUP::update()
00236 {
00237   static double u[2][9];
00238   int i;
00239   for (i = 0; i < nenu; i++) {
00240     const Vector &disp = theNodes[i]->getTrialDisp();
00241     u[0][i] = disp(0);
00242     u[1][i] = disp(1);
00243   }
00244   
00245   static Vector eps(3);
00246   
00247   int ret = 0;
00248   
00249   // Determine Jacobian for this integration point
00250   this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 
00251 
00252   // Loop over the integration points
00253   for (i = 0; i < nintu; i++) {
00254     
00255     // Interpolate strains
00256     //eps = B*u;
00257     //eps.addMatrixVector(0.0, B, u, 1.0);
00258     eps.Zero();
00259     for (int beta = 0; beta < nenu; beta++) {
00260       eps(0) += shgu[0][beta][i]*u[0][beta];
00261       eps(1) += shgu[1][beta][i]*u[1][beta];
00262       eps(2) += shgu[0][beta][i]*u[1][beta] + shgu[1][beta][i]*u[0][beta];
00263     }
00264     
00265     // Set the material strain
00266     ret += theMaterial[i]->setTrialStrain(eps);
00267   }
00268   
00269   return ret;
00270 }
00271 
00272 
00273 const Matrix&
00274 NineFourNodeQuadUP::getTangentStiff()
00275 {
00276   int i, j, j2, j2m1, ik, ib, jk, jb;
00277   static Matrix B(3,nenu*2);
00278   static Matrix BTDB(nenu*2,nenu*2);
00279 
00280   B.Zero(); 
00281   BTDB.Zero();
00282   K.Zero();
00283 
00284   // Determine Jacobian for this integration point
00285   this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 
00286   
00287   // Loop over the integration points
00288   for (i = 0; i < nintu; i++) {
00289     
00290     // Get the material tangent
00291     const Matrix &D = theMaterial[i]->getTangent();
00292     
00293         for (j=0; j<nenu; j++) {
00294                 j2 = j*2+1;
00295                 j2m1 = j*2;
00296         B(0,j2m1) = shgu[0][j][i];
00297                 B(0,j2)   = 0.;
00298                 B(1,j2m1) = 0.;
00299                 B(1,j2)   = shgu[1][j][i];
00300                 B(2,j2m1) = shgu[1][j][i];
00301                 B(2,j2)   = shgu[0][j][i];
00302     }
00303 
00304     // Perform numerical integration
00305     //K = K + (B^ D * B) * intWt(i)*intWt(j) * detJ;
00306     BTDB.addMatrixTripleProduct(1.0, B, D, dvolu[i]);
00307   }
00308 
00309   for (i = 0; i < nenu; i++) {
00310           if (i<nenp) ik = i*3;
00311       if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00312       ib = i*2;
00313 
00314           for (j = 0; j < nenu; j++) {
00315                   if (j<nenp) jk = j*3;
00316                   if (j>=nenp) jk = nenp*3 + (j-nenp)*2;
00317           jb = j*2;
00318 
00319           K(ik,jk) += BTDB(ib,jb);
00320                   K(ik+1,jk) += BTDB(ib+1,jb);
00321                   K(ik,jk+1) += BTDB(ib,jb+1);
00322                   K(ik+1,jk+1) += BTDB(ib+1,jb+1);
00323           }
00324   }
00325 
00326   return K;
00327 }
00328 
00329 
00330 const Matrix &NineFourNodeQuadUP::getInitialStiff () 
00331 {
00332   if (Ki != 0) return *Ki;
00333 
00334   int i, j, j2, j2m1, ik, ib, jk, jb;
00335   static Matrix B(3,nenu*2);
00336   static Matrix BTDB(nenu*2,nenu*2);
00337 
00338   B.Zero(); 
00339   BTDB.Zero();
00340   K.Zero();
00341 
00342   // Determine Jacobian for this integration point
00343   this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 
00344   
00345   // Loop over the integration points
00346   for (i = 0; i < nintu; i++) {
00347     
00348     // Get the material tangent
00349     const Matrix &D = theMaterial[i]->getInitialTangent();
00350     
00351         for (j=0; j<nenu; j++) {
00352                 j2 = j*2+1;
00353                 j2m1 = j*2;
00354         B(0,j2m1) = shgu[0][j][i];
00355                 B(0,j2)   = 0.;
00356                 B(1,j2m1) = 0.;
00357                 B(1,j2)   = shgu[1][j][i];
00358                 B(2,j2m1) = shgu[1][j][i];
00359                 B(2,j2)   = shgu[0][j][i];
00360     }
00361 
00362     // Perform numerical integration
00363     //K = K + (B^ D * B) * intWt(i)*intWt(j) * detJ;
00364     BTDB.addMatrixTripleProduct(1.0, B, D, dvolu[i]);
00365   }
00366 
00367   for (i = 0; i < nenu; i++) {
00368           if (i<nenp) ik = i*3;
00369       if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00370       ib = i*2;
00371 
00372           for (j = 0; j < nenu; j++) {
00373                   if (j<nenp) jk = j*3;
00374                   if (j>=nenp) jk = nenp*3 + (j-nenp)*2;
00375           jb = j*2;
00376           K(ik,jk) += BTDB(ib,jb);
00377                   K(ik+1,jk) += BTDB(ib+1,jb);
00378                   K(ik,jk+1) += BTDB(ib,jb+1);
00379                   K(ik+1,jk+1) += BTDB(ib+1,jb+1);
00380           }
00381   }
00382 
00383   Ki = new Matrix(K);
00384   if (Ki == 0) {
00385     opserr << "FATAL NineFourNodeQuadUP::getInitialStiff() -";
00386     opserr << "ran out of memory\n";
00387     exit(-1);
00388   }  
00389     
00390   return *Ki;
00391 }
00392 
00393 
00394 const Matrix&
00395 NineFourNodeQuadUP::getDamp()
00396 {
00397   static Matrix Kdamp(22,22);
00398   Kdamp.Zero();
00399   
00400   if (betaK != 0.0)
00401     Kdamp.addMatrix(1.0, this->getTangentStiff(), betaK);      
00402   if (betaK0 != 0.0)
00403     Kdamp.addMatrix(1.0, this->getInitialStiff(), betaK0);      
00404   if (betaKc != 0.0)
00405     Kdamp.addMatrix(1.0, *Kc, betaKc);      
00406 
00407   int i, j, m, ik, jk;
00408 
00409   if (alphaM != 0.0) {
00410          this->getMass();
00411 
00412      for (i = 0; i < nenu; i++) {
00413              if (i<nenp) ik = i*3;
00414          if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00415 
00416              for (j = 0; j < nenu; j++) {
00417                      if (j<nenp) jk = j*3;
00418                      if (j>=nenp) jk = nenp*3 + (j-nenp)*2;
00419 
00420              Kdamp(ik,jk) += K(ik,jk)*alphaM;
00421              Kdamp(ik+1,jk+1) += K(ik+1,jk+1)*alphaM;
00422                  }
00423      }  
00424   }
00425 
00426   // Determine Jacobian for this integration point
00427   this->globalShapeFunction(dvolq, wp, nintp, nenu, 2); 
00428   this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 
00429 
00430   // Compute coupling matrix
00431   for (i = 0; i < nenu; i++) {
00432            if (i<nenp) ik = i*3;
00433        if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00434 
00435        for (j = 0; j < nenp; j++) {
00436                jk = j*3+2;
00437 
00438            for (m = 0; m < nintp; m++) {
00439                    Kdamp(ik,jk) += -dvolq[m]*shgq[0][i][m]*shgp[2][j][m];
00440                    Kdamp(ik+1,jk) += -dvolq[m]*shgq[1][i][m]*shgp[2][j][m];
00441                    }
00442            Kdamp(jk,ik) = Kdamp(ik,jk);
00443            Kdamp(jk,ik+1) = Kdamp(ik+1,jk);
00444            }
00445   }
00446   
00447   // Compute permeability matrix
00448   for (i = 0; i < nenp; i++) {
00449        ik = i*3+2;
00450 
00451        for (j = 0; j < nenp; j++) {
00452                jk = j*3+2;
00453 
00454            for (m = 0; m < nintp; m++) {
00455                Kdamp(ik,jk) += - dvolp[m]*(perm[0]*shgp[0][i][m]*shgp[0][j][m] +
00456                                    perm[1]*shgp[1][i][m]*shgp[1][j][m]);
00457                    }
00458            }
00459   }
00460 
00461   K = Kdamp;
00462   return K;
00463 }
00464 
00465 const Matrix&
00466 NineFourNodeQuadUP::getMass()
00467 {
00468   K.Zero();
00469   
00470   int i, j, m, ik, jk;
00471   double Nrho;
00472   
00473   // Determine Jacobian for this integration point
00474   this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 
00475 
00476   // Compute consistent mass matrix
00477   for (i = 0; i < nenu; i++) {
00478           if (i<nenp) ik = i*3;
00479       if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00480 
00481           for (j = 0; j < nenu; j++) {
00482                   if (j<nenp) jk = j*3;
00483               if (j>=nenp) jk = nenp*3 + (j-nenp)*2;
00484 
00485           for (m = 0; m < nintu; m++) {
00486               Nrho = dvolu[m]*mixtureRho(m)*shgu[2][i][m]*shgu[2][j][m];
00487               K(ik,jk) += Nrho;
00488               K(ik+1,jk+1) += Nrho;
00489                   }
00490           }
00491   }
00492   
00493   // Compute compressibility matrix
00494   double oneOverKc = 1./kc;
00495   this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 
00496   
00497   for (i = 0; i < nenp; i++) {
00498        ik = i*3+2;
00499 
00500        for (j = 0; j < nenp; j++) {
00501                jk = j*3+2;
00502 
00503            for (m = 0; m < nintp; m++) {
00504                   K(ik,jk) += -dvolp[m]*oneOverKc*shgp[2][i][m]*shgp[2][j][m];
00505           }
00506     }
00507   }
00508 
00509   return K;
00510 }
00511 
00512 void
00513 NineFourNodeQuadUP::zeroLoad(void)
00514 {
00515         Q.Zero();
00516 
00517         return;
00518 }
00519 
00520 int 
00521 NineFourNodeQuadUP::addLoad(ElementalLoad *theLoad, double loadFactor)
00522 {
00523   opserr << "NineFourNodeQuadUP::addLoad - load type unknown for ele with tag: " << this->getTag() << "\n";
00524   return -1;
00525 }
00526 
00527 
00528 int 
00529 NineFourNodeQuadUP::addInertiaLoadToUnbalance(const Vector &accel)
00530 {
00531   // accel = uDotDotG (see EarthquakePattern.cpp)
00532   // Get R * accel from the nodes
00533 
00534   static Vector ra(22);
00535   int i, j, ik;
00536 
00537   ra.Zero();
00538 
00539   for (i=0; i<nenu; i++) {
00540       const Vector &Raccel = theNodes[i]->getRV(accel);
00541       if ((i<nenp && 3 != Raccel.Size()) || (i>=nenp && 2 != Raccel.Size())) {
00542          opserr << "NineFourNodeQuadUP::addInertiaLoadToUnbalance matrix and vector sizes are incompatable\n";
00543          return -1;
00544           }
00545 
00546           if (i<nenp) ik = i*3;
00547       if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00548       ra[ik] = Raccel(0);
00549       ra[ik+1] = Raccel(1);
00550   }
00551 
00552   // Compute mass matrix
00553   this->getMass();
00554   
00555   // Want to add ( - fact * M R * accel ) to unbalance 
00556   for (i = 0; i < 22; i++) {
00557     for (j = 0; j < 22; j++)
00558       Q(i) += -K(i,j)*ra[j];
00559   }
00560   
00561   return 0;
00562 }
00563 
00564 const Vector&
00565 NineFourNodeQuadUP::getResistingForce()
00566 {
00567   P.Zero();
00568   
00569   int i, j, jk;
00570 
00571   // Determine Jacobian for this integration point
00572   this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 
00573   this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 
00574   
00575   // Loop over the integration points
00576   for (i = 0; i < nintu; i++) {
00577 
00578     // Get material stress response
00579     const Vector &sigma = theMaterial[i]->getStress();
00580     
00581     // Perform numerical integration on internal force
00582     //P = P + (B^ sigma) * intWt(i)*intWt(j) * detJ;
00583     //P.addMatrixTransposeVector(1.0, B, sigma, intWt(i)*intWt(j)*detJ);
00584         for (j = 0; j < nenu; j++) {
00585                 if (j<nenp) jk = j*3;
00586             if (j>=nenp) jk = nenp*3 + (j-nenp)*2;
00587 
00588         P(jk) += dvolu[i]*(shgu[0][j][i]*sigma(0) + shgu[1][j][i]*sigma(2));      
00589         P(jk+1) += dvolu[i]*(shgu[1][j][i]*sigma(1) + shgu[0][j][i]*sigma(2));
00590       
00591         // Subtract equiv. body forces from the nodes
00592         //P = P - (N^ b) * intWt(i)*intWt(j) * detJ;
00593         //P.addMatrixTransposeVector(1.0, N, b, -intWt(i)*intWt(j)*detJ);
00594         double r = mixtureRho(i);
00595         P(jk) -= dvolu[i]*(shgu[2][j][i]*r*b[0]);
00596         P(jk+1) -= dvolu[i]*(shgu[2][j][i]*r*b[1]);
00597     }
00598   }
00599 //  opserr<<"K -body"<<P<<endln;
00600   
00601   // Subtract fluid body force
00602   for (j = 0; j < nenp; j++) {
00603          jk = j*3+2;
00604      for (i = 0; i < nintp; i++) {
00605         P(jk) += dvolp[i]*rho*(perm[0]*b[0]*shgp[0][j][i] +
00606                  perm[1]*b[1]*shgp[1][j][i]);
00607      }
00608   }
00609 
00610 //  opserr<<"K -fluid "<<P<<endln;
00611 
00612   // Subtract other external nodal loads ... P_res = P_int - P_ext
00613   //P = P - Q;
00614   P.addVector(1.0, Q, -1.0);
00615   
00616   return P;
00617 }
00618 
00619 const Vector&
00620 NineFourNodeQuadUP::getResistingForceIncInertia()
00621 {
00622   int i, j, ik;
00623   static double a[22];
00624 
00625   for (i=0; i<nenu; i++) {
00626     const Vector &accel = theNodes[i]->getTrialAccel();
00627     if ((i<nenp && 3 != accel.Size()) || (i>=nenp && 2 != accel.Size())) {
00628       opserr << "NineFourNodeQuadUP::getResistingForceIncInertia matrix and vector sizes are incompatable\n";
00629          return P;
00630     }
00631     
00632     if (i<nenp) ik = i*3;
00633     if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00634     a[ik] = accel(0);
00635       a[ik+1] = accel(1);
00636       if (i<nenp) a[ik+2] = accel(2);
00637   }
00638   
00639   // Compute the current resisting force
00640   this->getResistingForce();
00641   //  opserr<<"K "<<P<<endln;
00642   
00643   // Compute the mass matrix
00644   this->getMass();
00645   
00646   for (i = 0; i < 22; i++) {
00647     for (j = 0; j < 22; j++)
00648       P(i) += K(i,j)*a[j];
00649   }
00650 //  opserr<<"K+M "<<P<<endln; 
00651    
00652   for (i=0; i<nenu; i++) {
00653       const Vector &vel = theNodes[i]->getTrialVel();
00654       if ((i<nenp && 3 != vel.Size()) || (i>=nenp && 2 != vel.Size())) {
00655          opserr << "NineFourNodeQuadUP::getResistingForceIncInertia matrix and vector sizes are incompatable\n";
00656          return P;
00657       }
00658 
00659       if (i<nenp) ik = i*3;
00660       if (i>=nenp) ik = nenp*3 + (i-nenp)*2;
00661       a[ik] = vel(0);
00662       a[ik+1] = vel(1);
00663       if (i<nenp) a[ik+2] = vel(2);
00664   }
00665   
00666   this->getDamp();
00667   
00668   for (i = 0; i < 22; i++) {
00669     for (j = 0; j < 22; j++) {
00670       P(i) += K(i,j)*a[j];
00671     }
00672   }
00673 //  opserr<<"final "<<P<<endln;
00674   return P;
00675 }
00676 
00677 int
00678 NineFourNodeQuadUP::sendSelf(int commitTag, Channel &theChannel)
00679 {
00680   int res = 0;
00681   
00682   // note: we don't check for dataTag == 0 for Element
00683   // objects as that is taken care of in a commit by the Domain
00684   // object - don't want to have to do the check if sending data
00685   int dataTag = this->getDbTag();
00686   
00687   // Quad packs its data into a Vector and sends this to theChannel
00688         // along with its dbTag and the commitTag passed in the arguments
00689   static Vector data(10);
00690   data(0) = this->getTag();
00691   data(1) = thickness;
00692   data(2) = rho;
00693   data(3) = b[0];
00694   data(4) = b[1];
00695   //data(5) = pressure;
00696   data(6) = kc;
00697   data(7) = perm[0];
00698   data(8) = perm[1];
00699   data(9) = 9;
00700   
00701   
00702   res += theChannel.sendVector(dataTag, commitTag, data);
00703   if (res < 0) {
00704     opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send Vector\n";
00705     return res;
00706   }           
00707   
00708   // Quad then sends the tags of its four end nodes
00709   res += theChannel.sendID(dataTag, commitTag, connectedExternalNodes);
00710   if (res < 0) {
00711     opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send ID\n";
00712     return res;
00713   }
00714   
00715   // Now quad sends the ids of its materials
00716   int matDbTag;
00717   int numMats = 9;
00718   ID classTags(2*numMats);
00719   
00720   int i;
00721   for (i = 0; i < 9; i++) {
00722     classTags(i) = theMaterial[i]->getClassTag();
00723     matDbTag = theMaterial[i]->getDbTag();
00724     // NOTE: we do have to ensure that the material has a database
00725     // tag if we are sending to a database channel.
00726     if (matDbTag == 0) {
00727       matDbTag = theChannel.getDbTag();
00728       if (matDbTag != 0)
00729         theMaterial[i]->setDbTag(matDbTag);
00730     }
00731     classTags(i+numMats) = matDbTag;
00732   }
00733   
00734   res += theChannel.sendID(dataTag, commitTag, classTags);
00735   if (res < 0) {
00736     opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send ID\n";
00737     return res;
00738   }
00739   
00740   // Finally, quad asks its material objects to send themselves
00741   for (i = 0; i < 9; i++) {
00742     res += theMaterial[i]->sendSelf(commitTag, theChannel);
00743     if (res < 0) {
00744       opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send its Material\n";
00745       return res;
00746     }
00747   }
00748   
00749   return res;
00750 }
00751 
00752 int
00753 NineFourNodeQuadUP::recvSelf(int commitTag, Channel &theChannel,
00754                                                 FEM_ObjectBroker &theBroker)
00755 {
00756   int res = 0;
00757   
00758   int dataTag = this->getDbTag();
00759   
00760   // Quad creates a Vector, receives the Vector and then sets the 
00761   // internal data with the data in the Vector
00762   static Vector data(7);
00763   res += theChannel.recvVector(dataTag, commitTag, data);
00764   if (res < 0) {
00765     opserr << "WARNING NineFourNodeQuadUP::recvSelf() - failed to receive Vector\n";
00766     return res;
00767   }
00768   
00769   this->setTag((int)data(0));
00770   thickness = data(1);
00771   rho = data(2);
00772   b[0] = data(3);
00773   b[1] = data(4);
00774   //pressure = data(5);
00775   kc = data(6);
00776   perm[0] = data(7);
00777   perm[1] = data(8);
00778 
00779   // Quad now receives the tags of its four external nodes
00780   res += theChannel.recvID(dataTag, commitTag, connectedExternalNodes);
00781   if (res < 0) {
00782     opserr << "WARNING NineFourNodeQuadUP::recvSelf() - " << this->getTag() << " failed to receive ID\n";
00783     return res;
00784   }
00785 
00786   // Quad now receives the ids of its materials
00787   int newOrder = (int)data(9);
00788   int numMats = newOrder;
00789   ID classTags(2*numMats);
00790 
00791   res += theChannel.recvID(dataTag, commitTag, classTags);
00792   if (res < 0)  {
00793     opserr << "NineFourNodeQuadUP::recvSelf() - failed to recv ID data\n";
00794     return res;
00795   }    
00796 
00797   int i;
00798   
00799   // If the number of materials (quadrature order) is not the same,
00800   // delete the old materials, allocate new ones and then receive
00801   if (9 != newOrder) {
00802                 // Delete the materials
00803     for (i = 0; i < 9; i++) {
00804       if (theMaterial[i])
00805         delete theMaterial[i];
00806                 }
00807     if (theMaterial)
00808       delete [] theMaterial;
00809     
00810     // Allocate new materials
00811     theMaterial = new NDMaterial *[9];
00812     if (theMaterial == 0) {
00813       opserr << "NineFourNodeQuadUP::recvSelf() - Could not allocate NDMaterial* array\n";
00814       return -1;
00815     }
00816     for (i = 0; i < 9; i++) {
00817       int matClassTag = classTags(i);
00818       int matDbTag = classTags(i+numMats);
00819       // Allocate new material with the sent class tag
00820       theMaterial[i] = theBroker.getNewNDMaterial(matClassTag);
00821       if (theMaterial[i] == 0) {
00822         opserr << "NineFourNodeQuadUP::recvSelf() - Broker could not create NDMaterial of class type" << matClassTag << endln;
00823         return -1;
00824       }
00825       // Now receive materials into the newly allocated space
00826       theMaterial[i]->setDbTag(matDbTag);
00827       res += theMaterial[i]->recvSelf(commitTag, theChannel, theBroker);
00828       if (res < 0) {
00829         opserr << "NLBeamColumn3d::recvSelf() - material " << i << "failed to recv itself\n";
00830         return res;
00831       }
00832     }
00833   }
00834   // Number of materials is the same, receive materials into current space
00835   else {
00836     for (i = 0; i < 9; i++) {
00837       int matClassTag = classTags(i);
00838       int matDbTag = classTags(i+numMats);
00839       // Check that material is of the right type; if not,
00840       // delete it and create a new one of the right type
00841       if (theMaterial[i]->getClassTag() != matClassTag) {
00842         delete theMaterial[i];
00843         theMaterial[i] = theBroker.getNewNDMaterial(matClassTag);
00844         if (theMaterial[i] == 0) {
00845           opserr << "NineFourNodeQuadUP::recvSelf() - Broker could not create NDMaterial of class type " << matClassTag << endln;
00846           exit(-1);
00847         }
00848       }
00849       // Receive the material
00850       theMaterial[i]->setDbTag(matDbTag);
00851       res += theMaterial[i]->recvSelf(commitTag, theChannel, theBroker);
00852       if (res < 0) {
00853         opserr << "NineFourNodeQuadUP::recvSelf() - material " << i << "failed to recv itself\n";
00854         return res;
00855       }
00856     }
00857   }
00858   
00859   return res;
00860 }
00861 
00862 void
00863 NineFourNodeQuadUP::Print(OPS_Stream &s, int flag)
00864 {
00865         s << "\nNineFourNodeQuadUP, element id:  " << this->getTag() << endln;
00866         s << "\tConnected external nodes:  " << connectedExternalNodes;
00867         s << "\tthickness:  " << thickness << endln;
00868         s << "\tmass density:  " << rho << endln;
00869         //s << "\tsurface pressure:  " << pressure << endln;
00870         s << "\tbody forces:  " << b[0] << ' ' << b[1] << endln;
00871         theMaterial[0]->Print(s,flag);
00872         s << "\tStress (xx yy xy)" << endln;
00873         for (int i = 0; i < 9; i++)
00874                 s << "\t\tGauss point " << i+1 << ": " << theMaterial[i]->getStress();
00875 }
00876 
00877 int
00878 NineFourNodeQuadUP::displaySelf(Renderer &theViewer, int displayMode, float fact)
00879 {
00880     // first set the quantity to be displayed at the nodes;
00881     // if displayMode is 1 through 3 we will plot material stresses otherwise 0.0
00882 
00883     static Vector values(9);
00884 
00885     for (int j=0; j<9; j++)
00886            values(j) = 0.0;
00887 
00888     if (displayMode < 4 && displayMode > 0) {
00889         for (int i=0; i<9; i++) {
00890           const Vector &stress = theMaterial[i]->getStress();
00891           values(i) = stress(displayMode-1);
00892         }
00893     }
00894 
00895     // now  determine the end points of the quad based on
00896     // the display factor (a measure of the distorted image)
00897     // store this information in 4 3d vectors v1 through v4
00898     /*const Vector &end1Crd = nd1Ptr->getCrds();
00899     const Vector &end2Crd = nd2Ptr->getCrds();  
00900     const Vector &end3Crd = nd3Ptr->getCrds();  
00901     const Vector &end4Crd = nd4Ptr->getCrds();  
00902 
00903     const Vector &end1Disp = nd1Ptr->getDisp();
00904     const Vector &end2Disp = nd2Ptr->getDisp();
00905     const Vector &end3Disp = nd3Ptr->getDisp();
00906     const Vector &end4Disp = nd4Ptr->getDisp();
00907 
00908     static Matrix coords(4,3);
00909 
00910     for (int i = 0; i < 2; i++) {
00911       coords(0,i) = end1Crd(i) + end1Disp(i)*fact;
00912       coords(1,i) = end2Crd(i) + end2Disp(i)*fact;    
00913       coords(2,i) = end3Crd(i) + end3Disp(i)*fact;    
00914       coords(3,i) = end4Crd(i) + end4Disp(i)*fact;    
00915     }
00916 
00917     int error = 0;
00918 
00919     // finally we draw the element using drawPolygon
00920     error += theViewer.drawPolygon (coords, values);
00921 
00922     return error;*/
00923         return 0;
00924 }
00925 
00926 Response*
00927 NineFourNodeQuadUP::setResponse(const char **argv, int argc, Information &eleInfo)
00928 {
00929     if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0)
00930                 return new ElementResponse(this, 1, P);
00931     
00932     else if (strcmp(argv[0],"stiff") == 0 || strcmp(argv[0],"stiffness") == 0)
00933                 return new ElementResponse(this, 2, K);
00934 
00935     else if (strcmp(argv[0],"mass") == 0)
00936                 return new ElementResponse(this, 3, K);
00937 
00938     else if (strcmp(argv[0],"damp") == 0)
00939                 return new ElementResponse(this, 4, K);
00940 
00941         else if (strcmp(argv[0],"material") == 0 || strcmp(argv[0],"integrPoint") == 0) {
00942                 int pointNum = atoi(argv[1]);
00943                 if (pointNum > 0 && pointNum <= nenu)
00944                         return theMaterial[pointNum-1]->setResponse(&argv[2], argc-2, eleInfo);
00945             else 
00946                         return 0;
00947         }
00948  
00949     // otherwise response quantity is unknown for the quad class
00950     else
00951                 return 0;
00952 }
00953 
00954 int 
00955 NineFourNodeQuadUP::getResponse(int responseID, Information &eleInfo)
00956 {
00957         switch (responseID) {
00958       
00959                 case 1:
00960                         return eleInfo.setVector(this->getResistingForce());
00961       
00962                 case 2:
00963                         return eleInfo.setMatrix(this->getTangentStiff());
00964 
00965                 case 3:
00966                         return eleInfo.setMatrix(this->getMass());
00967 
00968                 case 4:
00969                         return eleInfo.setMatrix(this->getDamp());
00970 
00971                 default: 
00972                         return -1;
00973         }
00974 }
00975 
00976 int
00977 NineFourNodeQuadUP::setParameter(const char **argv, int argc, Information &info)
00978 {
00979         // quad mass density per unit volume
00980         if (strcmp(argv[0],"rho") == 0) {
00981                 info.theType = DoubleType;
00982                 info.theDouble = rho;
00983                 return 1;
00984         }
00985 
00986     // a material parameter
00987     else if (strcmp(argv[0],"material") == 0) {
00988                 int pointNum = atoi(argv[1]);
00989                 if (pointNum > 0 && pointNum <= nenu) {
00990                         int ok = theMaterial[pointNum-1]->setParameter(&argv[2], argc-2, info);
00991                         if (ok < 0)
00992                                 return -1;
00993                     else if (ok >= 0 && ok < 100)
00994                                 return pointNum*100 + ok;
00995                         else
00996                                 return -1;
00997                 }
00998             else 
00999                         return -1;
01000         }
01001     
01002     // otherwise parameter is unknown for the Truss class
01003     else
01004                 return -1;
01005 
01006 }
01007     
01008 int
01009 NineFourNodeQuadUP::updateParameter(int parameterID, Information &info)
01010 {
01011   switch (parameterID) {
01012     case -1:
01013       return -1;
01014       
01015         case 1:
01016                 rho = info.theDouble;
01017                 this->getMass();        // update mass matrix
01018                 return 0;
01019 
01020         default: 
01021                 if (parameterID >= 100) { // material parameter
01022                         int pointNum = parameterID/100;
01023                         if (pointNum > 0 && pointNum <= nenu)
01024                                 return theMaterial[pointNum-1]->updateParameter(parameterID-100*pointNum, info);
01025                         else
01026                                 return -1;
01027                 } else // unknown
01028                         return -1;
01029   }
01030 }
01031 
01032 void NineFourNodeQuadUP::globalShapeFunction(double *dvol, double *w, int nint, int nen, int mode)
01033 {
01034   static double coord[2][9], xs[2][2], det, temp;
01035   int i, j, k, m;
01036 
01037   for (i=0; i<3; i++) {
01038      for (j=0; j<nen; j++) {
01039         for (k=0; k<nint; k++) {
01040            if (mode==0) shgu[i][j][k] = shlu[i][j][k];
01041            if (mode==1) shgp[i][j][k] = shlp[i][j][k];
01042            if (mode==2) shgq[i][j][k] = shlq[i][j][k];
01043                 }
01044          }
01045   }
01046 
01047   for (i=0; i<nen; i++) {
01048          const Vector &coordRef = theNodes[i]->getCrds();
01049          coord[0][i] = coordRef(0);
01050          coord[1][i] = coordRef(1);
01051   }
01052 
01053   for (m=0; m<nint; m++) {
01054       for (i=0; i<2; i++) {
01055               for (j=0; j<2; j++) {
01056                    xs[i][j] = 0.0;
01057                    for (k=0; k<nen; k++) {
01058                     if (mode==0) xs[i][j] += coord[j][k]*shgu[i][k][m]; 
01059                     if (mode==1) xs[i][j] += coord[j][k]*shgp[i][k][m]; 
01060                     if (mode==2) xs[i][j] += coord[j][k]*shgq[i][k][m]; 
01061                }
01062                   }
01063           }
01064 
01065           det = xs[0][0]*xs[1][1] - xs[0][1]*xs[1][0];
01066 
01067           if (det < 0.0) {
01068           opserr << "WARNING NineFourNodeQuadUP: Determinant<=0 in tag " 
01069                      << this->getTag();
01070                   exit(-1);
01071           }
01072 
01073       for (i=0; i<nen; i++) {
01074           if (mode==0) { 
01075              temp = (shgu[0][i][m]*xs[1][1] - shgu[1][i][m]*xs[0][1])/det;
01076                      shgu[1][i][m] = (-shgu[0][i][m]*xs[1][0] + shgu[1][i][m]*xs[0][0])/det;
01077                      shgu[0][i][m] = temp;
01078           }
01079           if (mode==1) { 
01080              temp = (shgp[0][i][m]*xs[1][1] - shgp[1][i][m]*xs[0][1])/det;
01081                      shgp[1][i][m] = (-shgp[0][i][m]*xs[1][0] + shgp[1][i][m]*xs[0][0])/det;
01082                      shgp[0][i][m] = temp;
01083           }
01084           if (mode==2) { 
01085              temp = (shgq[0][i][m]*xs[1][1] - shgq[1][i][m]*xs[0][1])/det;
01086                      shgq[1][i][m] = (-shgq[0][i][m]*xs[1][0] + shgq[1][i][m]*xs[0][0])/det;
01087                      shgq[0][i][m] = temp;
01088           }
01089           }
01090 
01091           dvol[m] = w[m]*thickness*det;
01092 
01093   }  //end of m loop
01094 
01095 }
01096 
01097 
01098 void NineFourNodeQuadUP::shapeFunction(double *w, int nint, int nen, int mode)
01099 {
01100   static const double ra[] = {-0.5,0.5,0.5,-0.5,0.,0.5,0.,-0.5,0.};
01101   static const double sa[] = {-0.5,-0.5,0.5,0.5,-0.5,0.,0.5,0.,0.};
01102 
01103   double g, r, s, shl19, shl29, shl39, tempr, temps;
01104   int ic, il, is;
01105 
01106   g = 0.;
01107   if (nint == 4) {
01108           g=2./sqrt(3.0);
01109           w[0] = w[1] = w[2] = w[3] = 1.;
01110   }
01111   if (nint == 9) {
01112           g=2.*sqrt(3.0/5.0);
01113           w[0] = w[1] = w[2] = w[3] = 25./81.;
01114           w[4] = w[5] = w[6] = w[7] = 40./81.;
01115           w[8] = 64./81.;
01116   }
01117 
01118   for (int i=0; i<nint; i++) {
01119     r = g*ra[i];
01120         s = g*sa[i];
01121         shl19 = shl29 = shl39 = 0.;
01122         if (nen > 4) {
01123                 tempr = 1.-r*r;
01124                 temps = 1.-s*s;
01125                 if (nen == 9) {
01126                         if (mode==0) {
01127                            shlu[0][8][i] = -2.*r*temps;
01128                            shl19 = 0.5*shlu[0][8][i];
01129                            shlu[1][8][i] = -2.*s*tempr;
01130                            shl29 = 0.5*shlu[1][8][i];
01131                            shlu[2][8][i] = temps*tempr;
01132                            shl39 = 0.5*shlu[2][8][i];
01133                         }
01134                         if (mode==2) {
01135                            shlq[0][8][i] = -2.*r*temps;
01136                            shl19 = 0.5*shlq[0][8][i];
01137                            shlq[1][8][i] = -2.*s*tempr;
01138                            shl29 = 0.5*shlq[1][8][i];
01139                            shlq[2][8][i] = temps*tempr;
01140                            shl39 = 0.5*shlq[2][8][i];
01141                         }
01142         }
01143                 if (mode==0) {
01144                    shlu[0][4][i] = -r*(1.-s) - shl19;
01145                    shlu[1][4][i] = -0.5*tempr - shl29;
01146            shlu[2][4][i] = 0.5*tempr*(1.-s) - shl39;
01147                    shlu[0][5][i] = 0.5*temps - shl19;
01148                    shlu[1][5][i] = -s*(1.+r) - shl29;
01149            shlu[2][5][i] = 0.5*temps*(1.+r) - shl39;
01150                    shlu[0][6][i] = -r*(1.+s) - shl19;
01151                    shlu[1][6][i] = 0.5*tempr - shl29;
01152            shlu[2][6][i] = 0.5*tempr*(1.+s) - shl39;
01153                    shlu[0][7][i] = -0.5*temps - shl19;
01154                    shlu[1][7][i] = -s*(1.-r) - shl29;
01155            shlu[2][7][i] = 0.5*temps*(1.-r) - shl39;
01156                 }
01157                 if (mode==2) {
01158                    shlq[0][4][i] = -r*(1.-s) - shl19;
01159                    shlq[1][4][i] = -0.5*tempr - shl29;
01160            shlq[2][4][i] = 0.5*tempr*(1.-s) - shl39;
01161                    shlq[0][5][i] = 0.5*temps - shl19;
01162                    shlq[1][5][i] = -s*(1.+r) - shl29;
01163            shlq[2][5][i] = 0.5*temps*(1.+r) - shl39;
01164                    shlq[0][6][i] = -r*(1.+s) - shl19;
01165                    shlq[1][6][i] = 0.5*tempr - shl29;
01166            shlq[2][6][i] = 0.5*tempr*(1.+s) - shl39;
01167                    shlq[0][7][i] = -0.5*temps - shl19;
01168                    shlq[1][7][i] = -s*(1.-r) - shl29;
01169            shlq[2][7][i] = 0.5*temps*(1.-r) - shl39;
01170                 }
01171         }
01172 
01173     for (int k=0; k<4; k++) {
01174         tempr = 0.5 + ra[k]*r;
01175                 temps = 0.5 + sa[k]*s;
01176                 if (mode==0) {
01177                    shlu[0][k][i] = ra[k]*temps - 0.5*shl19;
01178                    shlu[1][k][i] = tempr*sa[k] - 0.5*shl29;
01179            shlu[2][k][i] = tempr*temps - 0.5*shl39;
01180                 }
01181                 if (mode==1) {
01182                    shlp[0][k][i] = ra[k]*temps - 0.5*shl19;
01183                    shlp[1][k][i] = tempr*sa[k] - 0.5*shl29;
01184            shlp[2][k][i] = tempr*temps - 0.5*shl39;
01185                 }
01186                 if (mode==2) {
01187                    shlq[0][k][i] = ra[k]*temps - 0.5*shl19;
01188                    shlq[1][k][i] = tempr*sa[k] - 0.5*shl29;
01189            shlq[2][k][i] = tempr*temps - 0.5*shl39;
01190                 }
01191         }
01192 
01193     if (nen > 4) {
01194         for (int m=4; m<8; m++) {
01195             ic = m - 4;
01196                         il = m - 3;
01197                         is = 1;
01198                         if (m==7) {
01199                 ic = 0;
01200                                 il = 3;
01201                                 is = 3;
01202                         }
01203             for (int j=ic; j<=il; j+=is) {
01204                         if (mode==0) {
01205                    shlu[0][j][i] = shlu[0][j][i] - 0.5*shlu[0][m][i];
01206                    shlu[1][j][i] = shlu[1][j][i] - 0.5*shlu[1][m][i];
01207                    shlu[2][j][i] = shlu[2][j][i] - 0.5*shlu[2][m][i];
01208                                 }
01209                         if (mode==2) {
01210                    shlq[0][j][i] = shlq[0][j][i] - 0.5*shlq[0][m][i];
01211                    shlq[1][j][i] = shlq[1][j][i] - 0.5*shlq[1][m][i];
01212                    shlq[2][j][i] = shlq[2][j][i] - 0.5*shlq[2][m][i];
01213                                 }
01214                         }
01215                 }  //end of m loop
01216         }
01217   }  //end of i loop
01218 }
01219 
01220 
01221 double NineFourNodeQuadUP::mixtureRho(int i)
01222 {
01223   double rhoi, e, n;
01224         rhoi= theMaterial[i]->getRho();
01225         e = 0.7;  //theMaterial[i]->getVoidRatio();
01226   n = e / (1.0 + e);
01227   //return n * rho + (1.0-n) * rhoi;
01228         return rhoi;
01229 }
01230

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