OpenSees » SRC » element » UP-ucsd

FourNodeQuadUP.cpp

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

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