Nine_Four_Node_QuadUP.cppGo to the documentation of this file.00001 00002 00003 // Description: This file contains the class definition for // 00004 00005 // NineFourNodeQuadUP, a 9-4-node (9 node for solid and 4 node for fluid) // 00006 00007 // plane strain element for solid-fluid fully coupled analysis. This // 00008 00009 // implementation is a simplified u-p formulation of Biot theory // 00010 00011 // (u - solid displacement, p - fluid pressure). Each element node has two // 00012 00013 // DOFs for u and 1 DOF for p. // 00014 00015 // // 00016 00017 // Written by Zhaohui Yang (March 2004) // 00018 00019 // // 00020 00022 00023 00024 00025 // $Revision: 1.4 $ 00026 00027 // $Date: 2006/09/05 21:04:59 $ 00028 00029 // $Source: /usr/local/cvs/OpenSees/SRC/element/UP-ucsd/Nine_Four_Node_QuadUP.cpp,v $ 00030 00031 00032 00033 #include <Nine_Four_Node_QuadUP.h> 00034 00035 #include <Node.h> 00036 00037 #include <NDMaterial.h> 00038 00039 #include <Matrix.h> 00040 00041 #include <Vector.h> 00042 00043 #include <ID.h> 00044 00045 #include <Renderer.h> 00046 00047 #include <Domain.h> 00048 00049 #include <string.h> 00050 00051 #include <Information.h> 00052 #include <Parameter.h> 00053 00054 #include <Channel.h> 00055 00056 #include <FEM_ObjectBroker.h> 00057 00058 #include <ElementResponse.h> 00059 00060 00061 00062 Matrix NineFourNodeQuadUP::K(22,22); 00063 00064 Vector NineFourNodeQuadUP::P(22); 00065 00066 double NineFourNodeQuadUP::shgu[3][9][9]; 00067 00068 double NineFourNodeQuadUP::shgp[3][4][4]; 00069 00070 double NineFourNodeQuadUP::shgq[3][9][4]; 00071 00072 double NineFourNodeQuadUP::shlu[3][9][9]; 00073 00074 double NineFourNodeQuadUP::shlp[3][4][4]; 00075 00076 double NineFourNodeQuadUP::shlq[3][9][4]; 00077 00078 double NineFourNodeQuadUP::wu[9]; 00079 00080 double NineFourNodeQuadUP::wp[4]; 00081 00082 double NineFourNodeQuadUP::dvolu[9]; 00083 00084 double NineFourNodeQuadUP::dvolp[4]; 00085 00086 double NineFourNodeQuadUP::dvolq[4]; 00087 00088 const int NineFourNodeQuadUP::nintu=9; 00089 00090 const int NineFourNodeQuadUP::nintp=4; 00091 00092 const int NineFourNodeQuadUP::nenu=9; 00093 00094 const int NineFourNodeQuadUP::nenp=4; 00095 00096 00097 00098 NineFourNodeQuadUP::NineFourNodeQuadUP(int tag, 00099 00100 int nd1, int nd2, int nd3, int nd4,int nd5, int nd6, int nd7, int nd8,int nd9, 00101 00102 NDMaterial &m, const char *type, double t, double bulk, double r, 00103 00104 double p1, double p2, double b1, double b2) 00105 00106 :Element (tag, ELE_TAG_Nine_Four_Node_QuadUP), 00107 00108 theMaterial(0), connectedExternalNodes(9), 00109 00110 Ki(0), Q(22), thickness(t), kc(bulk), rho(r) 00111 00112 { 00113 00114 this->shapeFunction(wu, nintu, nenu, 0); 00115 00116 /* for( int L = 0; L < nintu; L++) { 00117 00118 for( int j = 0; j < nenu; j++) { 00119 00120 printf("%5d %5d %15.6e %15.6e %15.6e\n", L+1, j+1, 00121 00122 shlu[0][j][L],shlu[1][j][L],shlu[2][j][L]); 00123 00124 } 00125 00126 } 00127 00128 exit(-1); 00129 00130 */ 00131 00132 this->shapeFunction(wp, nintp, nenp, 1); 00133 00134 this->shapeFunction(wp, nintp, nenu, 2); 00135 00136 00137 00138 // Body forces 00139 00140 b[0] = b1; 00141 00142 b[1] = b2; 00143 00144 // Permeabilities 00145 00146 perm[0] = p1; 00147 00148 perm[1] = p2; 00149 00150 00151 00152 // Allocate arrays of pointers to NDMaterials 00153 00154 theMaterial = new NDMaterial *[nintu]; 00155 00156 00157 00158 if (theMaterial == 0) { 00159 00160 opserr << "NineFourNodeQuadUP::NineFourNodeQuadUP - failed allocate material model pointer\n"; 00161 00162 exit(-1); 00163 00164 } 00165 00166 00167 00168 for (int i = 0; i < nintu; i++) { 00169 00170 00171 00172 // Get copies of the material model for each integration point 00173 00174 theMaterial[i] = m.getCopy(type); 00175 00176 00177 00178 // Check allocation 00179 00180 if (theMaterial[i] == 0) { 00181 00182 opserr << "NineFourNodeQuadUP::NineFourNodeQuadUP -- failed to get a copy of material model\n"; 00183 00184 exit(-1); 00185 00186 } 00187 00188 } 00189 00190 00191 00192 // Set connected external node IDs 00193 00194 connectedExternalNodes(0) = nd1; 00195 00196 connectedExternalNodes(1) = nd2; 00197 00198 connectedExternalNodes(2) = nd3; 00199 00200 connectedExternalNodes(3) = nd4; 00201 00202 connectedExternalNodes(4) = nd5; 00203 00204 connectedExternalNodes(5) = nd6; 00205 00206 connectedExternalNodes(6) = nd7; 00207 00208 connectedExternalNodes(7) = nd8; 00209 00210 connectedExternalNodes(8) = nd9; 00211 00212 } 00213 00214 00215 00216 00217 00218 NineFourNodeQuadUP::NineFourNodeQuadUP() 00219 00220 :Element (0,ELE_TAG_Nine_Four_Node_QuadUP), 00221 00222 theMaterial(0), connectedExternalNodes(9), 00223 00224 Ki(0), Q(22), thickness(0.0), kc(0.0), rho(0.0) 00225 00226 { 00227 00228 this->shapeFunction(wu, nintu, nenu, 0); 00229 00230 this->shapeFunction(wp, nintp, nenp, 1); 00231 00232 this->shapeFunction(wp, nintp, nenu, 2); 00233 00234 } 00235 00236 00237 00238 NineFourNodeQuadUP::~NineFourNodeQuadUP() 00239 00240 { 00241 00242 for (int i = 0; i < nintu; i++) { 00243 00244 if (theMaterial[i]) 00245 00246 delete theMaterial[i]; 00247 00248 } 00249 00250 00251 00252 // Delete the array of pointers to NDMaterial pointer arrays 00253 00254 if (theMaterial) 00255 00256 delete [] theMaterial; 00257 00258 00259 00260 for (int ii = 0; ii < nenu; ii++) theNodes[ii] = 0 ; 00261 00262 00263 00264 if (Ki != 0) 00265 00266 delete Ki; 00267 00268 } 00269 00270 00271 00272 int 00273 00274 NineFourNodeQuadUP::getNumExternalNodes() const 00275 00276 { 00277 00278 return nenu; 00279 00280 } 00281 00282 00283 00284 const ID& 00285 00286 NineFourNodeQuadUP::getExternalNodes() 00287 00288 { 00289 00290 return connectedExternalNodes; 00291 00292 } 00293 00294 00295 00296 Node ** 00297 00298 NineFourNodeQuadUP::getNodePtrs() 00299 00300 { 00301 00302 return theNodes; 00303 00304 } 00305 00306 00307 00308 int 00309 00310 NineFourNodeQuadUP::getNumDOF() 00311 00312 { 00313 00314 return 22; 00315 00316 } 00317 00318 00319 00320 void 00321 00322 NineFourNodeQuadUP::setDomain(Domain *theDomain) 00323 00324 { 00325 00326 // Check Domain is not null - invoked when object removed from a domain 00327 00328 if (theDomain == 0) { 00329 00330 for (int i=0; i<nenu; i++) theNodes[i] = 0; 00331 00332 return; 00333 00334 } 00335 00336 00337 00338 int i; 00339 00340 for (i=0; i<nenu; i++) { 00341 00342 theNodes[i] = theDomain->getNode(connectedExternalNodes(i)); 00343 00344 if (theNodes[i] == 0) { 00345 00346 opserr << "FATAL ERROR NineFourNodeQuadUP, node not found in domain, tag " 00347 00348 << this->getTag(); 00349 00350 return; 00351 00352 } 00353 00354 } 00355 00356 00357 00358 int dof; 00359 00360 for (i=0; i<nenu; i++) { 00361 00362 dof = theNodes[i]->getNumberDOF(); 00363 00364 if ((i<nenp && dof != 3) || (i>=nenp && dof != 2)) { 00365 00366 opserr << "FATAL ERROR NineFourNodeQuadUP, has wrong number of DOFs at its nodes " 00367 00368 << this->getTag(); 00369 00370 return; 00371 00372 } 00373 00374 } 00375 00376 00377 00378 this->DomainComponent::setDomain(theDomain); 00379 00380 } 00381 00382 00383 00384 00385 00386 int 00387 00388 NineFourNodeQuadUP::commitState() 00389 00390 { 00391 00392 int retVal = 0; 00393 00394 00395 00396 // call element commitState to do any base class stuff 00397 00398 if ((retVal = this->Element::commitState()) != 0) { 00399 00400 opserr << "Nine_Four_Node_Quad_UP::commitState () - failed in base class"; 00401 00402 } 00403 00404 00405 00406 // Loop over the integration points and commit the material states 00407 00408 for (int i = 0; i < nintu; i++) 00409 00410 retVal += theMaterial[i]->commitState(); 00411 00412 00413 00414 return retVal; 00415 00416 } 00417 00418 00419 00420 int 00421 00422 NineFourNodeQuadUP::revertToLastCommit() 00423 00424 { 00425 00426 int retVal = 0; 00427 00428 00429 00430 // Loop over the integration points and revert to last committed state 00431 00432 for (int i = 0; i < nintu; i++) 00433 00434 retVal += theMaterial[i]->revertToLastCommit(); 00435 00436 00437 00438 return retVal; 00439 00440 } 00441 00442 00443 00444 int 00445 00446 NineFourNodeQuadUP::revertToStart() 00447 00448 { 00449 00450 int retVal = 0; 00451 00452 00453 00454 // Loop over the integration points and revert states to start 00455 00456 for (int i = 0; i < nintu; i++) 00457 00458 retVal += theMaterial[i]->revertToStart(); 00459 00460 00461 00462 return retVal; 00463 00464 } 00465 00466 00467 00468 int 00469 00470 NineFourNodeQuadUP::update() 00471 00472 { 00473 00474 static double u[2][9]; 00475 00476 int i; 00477 00478 for (i = 0; i < nenu; i++) { 00479 00480 const Vector &disp = theNodes[i]->getTrialDisp(); 00481 00482 u[0][i] = disp(0); 00483 00484 u[1][i] = disp(1); 00485 00486 } 00487 00488 00489 00490 static Vector eps(3); 00491 00492 00493 00494 int ret = 0; 00495 00496 00497 00498 // Determine Jacobian for this integration point 00499 00500 this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 00501 00502 00503 00504 // Loop over the integration points 00505 00506 for (i = 0; i < nintu; i++) { 00507 00508 00509 00510 // Interpolate strains 00511 00512 //eps = B*u; 00513 00514 //eps.addMatrixVector(0.0, B, u, 1.0); 00515 00516 eps.Zero(); 00517 00518 for (int beta = 0; beta < nenu; beta++) { 00519 00520 eps(0) += shgu[0][beta][i]*u[0][beta]; 00521 00522 eps(1) += shgu[1][beta][i]*u[1][beta]; 00523 00524 eps(2) += shgu[0][beta][i]*u[1][beta] + shgu[1][beta][i]*u[0][beta]; 00525 00526 } 00527 00528 00529 00530 // Set the material strain 00531 00532 ret += theMaterial[i]->setTrialStrain(eps); 00533 00534 } 00535 00536 00537 00538 return ret; 00539 00540 } 00541 00542 00543 00544 00545 00546 const Matrix& 00547 00548 NineFourNodeQuadUP::getTangentStiff() 00549 00550 { 00551 00552 int i, j, j2, j2m1, ik, ib, jk, jb; 00553 00554 static Matrix B(3,nenu*2); 00555 00556 static Matrix BTDB(nenu*2,nenu*2); 00557 00558 00559 00560 B.Zero(); 00561 00562 BTDB.Zero(); 00563 00564 K.Zero(); 00565 00566 00567 00568 // Determine Jacobian for this integration point 00569 00570 this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 00571 00572 00573 00574 // Loop over the integration points 00575 00576 for (i = 0; i < nintu; i++) { 00577 00578 00579 00580 // Get the material tangent 00581 00582 const Matrix &D = theMaterial[i]->getTangent(); 00583 00584 00585 00586 for (j=0; j<nenu; j++) { 00587 00588 j2 = j*2+1; 00589 00590 j2m1 = j*2; 00591 00592 B(0,j2m1) = shgu[0][j][i]; 00593 00594 B(0,j2) = 0.; 00595 00596 B(1,j2m1) = 0.; 00597 00598 B(1,j2) = shgu[1][j][i]; 00599 00600 B(2,j2m1) = shgu[1][j][i]; 00601 00602 B(2,j2) = shgu[0][j][i]; 00603 00604 } 00605 00606 00607 00608 // Perform numerical integration 00609 00610 //K = K + (B^ D * B) * intWt(i)*intWt(j) * detJ; 00611 00612 BTDB.addMatrixTripleProduct(1.0, B, D, dvolu[i]); 00613 00614 } 00615 00616 00617 00618 for (i = 0; i < nenu; i++) { 00619 00620 if (i<nenp) ik = i*3; 00621 00622 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 00623 00624 ib = i*2; 00625 00626 00627 00628 for (j = 0; j < nenu; j++) { 00629 00630 if (j<nenp) jk = j*3; 00631 00632 if (j>=nenp) jk = nenp*3 + (j-nenp)*2; 00633 00634 jb = j*2; 00635 00636 00637 00638 K(ik,jk) += BTDB(ib,jb); 00639 00640 K(ik+1,jk) += BTDB(ib+1,jb); 00641 00642 K(ik,jk+1) += BTDB(ib,jb+1); 00643 00644 K(ik+1,jk+1) += BTDB(ib+1,jb+1); 00645 00646 } 00647 00648 } 00649 00650 00651 00652 return K; 00653 00654 } 00655 00656 00657 00658 00659 00660 const Matrix &NineFourNodeQuadUP::getInitialStiff () 00661 00662 { 00663 00664 if (Ki != 0) return *Ki; 00665 00666 00667 00668 int i, j, j2, j2m1, ik, ib, jk, jb; 00669 00670 static Matrix B(3,nenu*2); 00671 00672 static Matrix BTDB(nenu*2,nenu*2); 00673 00674 00675 00676 B.Zero(); 00677 00678 BTDB.Zero(); 00679 00680 K.Zero(); 00681 00682 00683 00684 // Determine Jacobian for this integration point 00685 00686 this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 00687 00688 00689 00690 // Loop over the integration points 00691 00692 for (i = 0; i < nintu; i++) { 00693 00694 00695 00696 // Get the material tangent 00697 00698 const Matrix &D = theMaterial[i]->getInitialTangent(); 00699 00700 00701 00702 for (j=0; j<nenu; j++) { 00703 00704 j2 = j*2+1; 00705 00706 j2m1 = j*2; 00707 00708 B(0,j2m1) = shgu[0][j][i]; 00709 00710 B(0,j2) = 0.; 00711 00712 B(1,j2m1) = 0.; 00713 00714 B(1,j2) = shgu[1][j][i]; 00715 00716 B(2,j2m1) = shgu[1][j][i]; 00717 00718 B(2,j2) = shgu[0][j][i]; 00719 00720 } 00721 00722 00723 00724 // Perform numerical integration 00725 00726 //K = K + (B^ D * B) * intWt(i)*intWt(j) * detJ; 00727 00728 BTDB.addMatrixTripleProduct(1.0, B, D, dvolu[i]); 00729 00730 } 00731 00732 00733 00734 for (i = 0; i < nenu; i++) { 00735 00736 if (i<nenp) ik = i*3; 00737 00738 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 00739 00740 ib = i*2; 00741 00742 00743 00744 for (j = 0; j < nenu; j++) { 00745 00746 if (j<nenp) jk = j*3; 00747 00748 if (j>=nenp) jk = nenp*3 + (j-nenp)*2; 00749 00750 jb = j*2; 00751 00752 K(ik,jk) += BTDB(ib,jb); 00753 00754 K(ik+1,jk) += BTDB(ib+1,jb); 00755 00756 K(ik,jk+1) += BTDB(ib,jb+1); 00757 00758 K(ik+1,jk+1) += BTDB(ib+1,jb+1); 00759 00760 } 00761 00762 } 00763 00764 00765 00766 Ki = new Matrix(K); 00767 00768 if (Ki == 0) { 00769 00770 opserr << "FATAL NineFourNodeQuadUP::getInitialStiff() -"; 00771 00772 opserr << "ran out of memory\n"; 00773 00774 exit(-1); 00775 00776 } 00777 00778 00779 00780 return *Ki; 00781 00782 } 00783 00784 00785 00786 00787 00788 const Matrix& 00789 00790 NineFourNodeQuadUP::getDamp() 00791 00792 { 00793 00794 static Matrix Kdamp(22,22); 00795 00796 Kdamp.Zero(); 00797 00798 00799 00800 if (betaK != 0.0) 00801 00802 Kdamp.addMatrix(1.0, this->getTangentStiff(), betaK); 00803 00804 if (betaK0 != 0.0) 00805 00806 Kdamp.addMatrix(1.0, this->getInitialStiff(), betaK0); 00807 00808 if (betaKc != 0.0) 00809 00810 Kdamp.addMatrix(1.0, *Kc, betaKc); 00811 00812 00813 00814 int i, j, m, ik, jk; 00815 00816 00817 00818 if (alphaM != 0.0) { 00819 00820 this->getMass(); 00821 00822 00823 00824 for (i = 0; i < nenu; i++) { 00825 00826 if (i<nenp) ik = i*3; 00827 00828 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 00829 00830 00831 00832 for (j = 0; j < nenu; j++) { 00833 00834 if (j<nenp) jk = j*3; 00835 00836 if (j>=nenp) jk = nenp*3 + (j-nenp)*2; 00837 00838 00839 00840 Kdamp(ik,jk) += K(ik,jk)*alphaM; 00841 00842 Kdamp(ik+1,jk+1) += K(ik+1,jk+1)*alphaM; 00843 00844 } 00845 00846 } 00847 00848 } 00849 00850 00851 00852 // Determine Jacobian for this integration point 00853 00854 this->globalShapeFunction(dvolq, wp, nintp, nenu, 2); 00855 00856 this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 00857 00858 00859 00860 // Compute coupling matrix 00861 00862 for (i = 0; i < nenu; i++) { 00863 00864 if (i<nenp) ik = i*3; 00865 00866 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 00867 00868 00869 00870 for (j = 0; j < nenp; j++) { 00871 00872 jk = j*3+2; 00873 00874 00875 00876 for (m = 0; m < nintp; m++) { 00877 00878 Kdamp(ik,jk) += -dvolq[m]*shgq[0][i][m]*shgp[2][j][m]; 00879 00880 Kdamp(ik+1,jk) += -dvolq[m]*shgq[1][i][m]*shgp[2][j][m]; 00881 00882 } 00883 00884 Kdamp(jk,ik) = Kdamp(ik,jk); 00885 00886 Kdamp(jk,ik+1) = Kdamp(ik+1,jk); 00887 00888 } 00889 00890 } 00891 00892 00893 00894 // Compute permeability matrix 00895 00896 for (i = 0; i < nenp; i++) { 00897 00898 ik = i*3+2; 00899 00900 00901 00902 for (j = 0; j < nenp; j++) { 00903 00904 jk = j*3+2; 00905 00906 00907 00908 for (m = 0; m < nintp; m++) { 00909 00910 Kdamp(ik,jk) += - dvolp[m]*(perm[0]*shgp[0][i][m]*shgp[0][j][m] + 00911 00912 perm[1]*shgp[1][i][m]*shgp[1][j][m]); 00913 00914 } 00915 00916 } 00917 00918 } 00919 00920 00921 00922 K = Kdamp; 00923 00924 return K; 00925 00926 } 00927 00928 00929 00930 const Matrix& 00931 00932 NineFourNodeQuadUP::getMass() 00933 00934 { 00935 00936 K.Zero(); 00937 00938 00939 00940 int i, j, m, ik, jk; 00941 00942 double Nrho; 00943 00944 00945 00946 // Determine Jacobian for this integration point 00947 00948 this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 00949 00950 00951 00952 // Compute consistent mass matrix 00953 00954 for (i = 0; i < nenu; i++) { 00955 00956 if (i<nenp) ik = i*3; 00957 00958 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 00959 00960 00961 00962 for (j = 0; j < nenu; j++) { 00963 00964 if (j<nenp) jk = j*3; 00965 00966 if (j>=nenp) jk = nenp*3 + (j-nenp)*2; 00967 00968 00969 00970 for (m = 0; m < nintu; m++) { 00971 00972 Nrho = dvolu[m]*mixtureRho(m)*shgu[2][i][m]*shgu[2][j][m]; 00973 00974 K(ik,jk) += Nrho; 00975 00976 K(ik+1,jk+1) += Nrho; 00977 00978 } 00979 00980 } 00981 00982 } 00983 00984 00985 00986 // Compute compressibility matrix 00987 00988 double oneOverKc = 1./kc; 00989 00990 this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 00991 00992 00993 00994 for (i = 0; i < nenp; i++) { 00995 00996 ik = i*3+2; 00997 00998 00999 01000 for (j = 0; j < nenp; j++) { 01001 01002 jk = j*3+2; 01003 01004 01005 01006 for (m = 0; m < nintp; m++) { 01007 01008 K(ik,jk) += -dvolp[m]*oneOverKc*shgp[2][i][m]*shgp[2][j][m]; 01009 01010 } 01011 01012 } 01013 01014 } 01015 01016 01017 01018 return K; 01019 01020 } 01021 01022 01023 01024 void 01025 01026 NineFourNodeQuadUP::zeroLoad(void) 01027 01028 { 01029 01030 Q.Zero(); 01031 01032 01033 01034 return; 01035 01036 } 01037 01038 01039 01040 int 01041 01042 NineFourNodeQuadUP::addLoad(ElementalLoad *theLoad, double loadFactor) 01043 01044 { 01045 01046 opserr << "NineFourNodeQuadUP::addLoad - load type unknown for ele with tag: " << this->getTag() << "\n"; 01047 01048 return -1; 01049 01050 } 01051 01052 01053 01054 01055 01056 int 01057 01058 NineFourNodeQuadUP::addInertiaLoadToUnbalance(const Vector &accel) 01059 01060 { 01061 01062 // accel = uDotDotG (see EarthquakePattern.cpp) 01063 01064 // Get R * accel from the nodes 01065 01066 01067 01068 static Vector ra(22); 01069 01070 int i, j, ik; 01071 01072 01073 01074 ra.Zero(); 01075 01076 01077 01078 for (i=0; i<nenu; i++) { 01079 01080 const Vector &Raccel = theNodes[i]->getRV(accel); 01081 01082 if ((i<nenp && 3 != Raccel.Size()) || (i>=nenp && 2 != Raccel.Size())) { 01083 01084 opserr << "NineFourNodeQuadUP::addInertiaLoadToUnbalance matrix and vector sizes are incompatable\n"; 01085 01086 return -1; 01087 01088 } 01089 01090 01091 01092 if (i<nenp) ik = i*3; 01093 01094 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 01095 01096 ra[ik] = Raccel(0); 01097 01098 ra[ik+1] = Raccel(1); 01099 01100 } 01101 01102 01103 01104 // Compute mass matrix 01105 01106 this->getMass(); 01107 01108 01109 01110 // Want to add ( - fact * M R * accel ) to unbalance 01111 01112 for (i = 0; i < 22; i++) { 01113 01114 for (j = 0; j < 22; j++) 01115 01116 Q(i) += -K(i,j)*ra[j]; 01117 01118 } 01119 01120 01121 01122 return 0; 01123 01124 } 01125 01126 01127 01128 const Vector& 01129 01130 NineFourNodeQuadUP::getResistingForce() 01131 01132 { 01133 01134 P.Zero(); 01135 01136 01137 01138 int i, j, jk; 01139 01140 01141 01142 // Determine Jacobian for this integration point 01143 01144 this->globalShapeFunction(dvolu, wu, nintu, nenu, 0); 01145 01146 this->globalShapeFunction(dvolp, wp, nintp, nenp, 1); 01147 01148 01149 01150 // Loop over the integration points 01151 01152 for (i = 0; i < nintu; i++) { 01153 01154 01155 01156 // Get material stress response 01157 01158 const Vector &sigma = theMaterial[i]->getStress(); 01159 01160 01161 01162 // Perform numerical integration on internal force 01163 01164 //P = P + (B^ sigma) * intWt(i)*intWt(j) * detJ; 01165 01166 //P.addMatrixTransposeVector(1.0, B, sigma, intWt(i)*intWt(j)*detJ); 01167 01168 for (j = 0; j < nenu; j++) { 01169 01170 if (j<nenp) jk = j*3; 01171 01172 if (j>=nenp) jk = nenp*3 + (j-nenp)*2; 01173 01174 01175 01176 P(jk) += dvolu[i]*(shgu[0][j][i]*sigma(0) + shgu[1][j][i]*sigma(2)); 01177 01178 P(jk+1) += dvolu[i]*(shgu[1][j][i]*sigma(1) + shgu[0][j][i]*sigma(2)); 01179 01180 01181 01182 // Subtract equiv. body forces from the nodes 01183 01184 //P = P - (N^ b) * intWt(i)*intWt(j) * detJ; 01185 01186 //P.addMatrixTransposeVector(1.0, N, b, -intWt(i)*intWt(j)*detJ); 01187 01188 double r = mixtureRho(i); 01189 01190 P(jk) -= dvolu[i]*(shgu[2][j][i]*r*b[0]); 01191 01192 P(jk+1) -= dvolu[i]*(shgu[2][j][i]*r*b[1]); 01193 01194 } 01195 01196 } 01197 01198 // opserr<<"K -body"<<P<<endln; 01199 01200 01201 01202 // Subtract fluid body force 01203 01204 for (j = 0; j < nenp; j++) { 01205 01206 jk = j*3+2; 01207 01208 for (i = 0; i < nintp; i++) { 01209 01210 P(jk) += dvolp[i]*rho*(perm[0]*b[0]*shgp[0][j][i] + 01211 01212 perm[1]*b[1]*shgp[1][j][i]); 01213 01214 } 01215 01216 } 01217 01218 01219 01220 // opserr<<"K -fluid "<<P<<endln; 01221 01222 01223 01224 // Subtract other external nodal loads ... P_res = P_int - P_ext 01225 01226 //P = P - Q; 01227 01228 P.addVector(1.0, Q, -1.0); 01229 01230 01231 01232 return P; 01233 01234 } 01235 01236 01237 01238 const Vector& 01239 01240 NineFourNodeQuadUP::getResistingForceIncInertia() 01241 01242 { 01243 01244 int i, j, ik; 01245 01246 static double a[22]; 01247 01248 01249 01250 for (i=0; i<nenu; i++) { 01251 01252 const Vector &accel = theNodes[i]->getTrialAccel(); 01253 01254 if ((i<nenp && 3 != accel.Size()) || (i>=nenp && 2 != accel.Size())) { 01255 01256 opserr << "NineFourNodeQuadUP::getResistingForceIncInertia matrix and vector sizes are incompatable\n"; 01257 01258 return P; 01259 01260 } 01261 01262 01263 01264 if (i<nenp) ik = i*3; 01265 01266 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 01267 01268 a[ik] = accel(0); 01269 01270 a[ik+1] = accel(1); 01271 01272 if (i<nenp) a[ik+2] = accel(2); 01273 01274 } 01275 01276 01277 01278 // Compute the current resisting force 01279 01280 this->getResistingForce(); 01281 01282 // opserr<<"K "<<P<<endln; 01283 01284 01285 01286 // Compute the mass matrix 01287 01288 this->getMass(); 01289 01290 01291 01292 for (i = 0; i < 22; i++) { 01293 01294 for (j = 0; j < 22; j++) 01295 01296 P(i) += K(i,j)*a[j]; 01297 01298 } 01299 01300 // opserr<<"K+M "<<P<<endln; 01301 01302 01303 01304 for (i=0; i<nenu; i++) { 01305 01306 const Vector &vel = theNodes[i]->getTrialVel(); 01307 01308 if ((i<nenp && 3 != vel.Size()) || (i>=nenp && 2 != vel.Size())) { 01309 01310 opserr << "NineFourNodeQuadUP::getResistingForceIncInertia matrix and vector sizes are incompatable\n"; 01311 01312 return P; 01313 01314 } 01315 01316 01317 01318 if (i<nenp) ik = i*3; 01319 01320 if (i>=nenp) ik = nenp*3 + (i-nenp)*2; 01321 01322 a[ik] = vel(0); 01323 01324 a[ik+1] = vel(1); 01325 01326 if (i<nenp) a[ik+2] = vel(2); 01327 01328 } 01329 01330 01331 01332 this->getDamp(); 01333 01334 01335 01336 for (i = 0; i < 22; i++) { 01337 01338 for (j = 0; j < 22; j++) { 01339 01340 P(i) += K(i,j)*a[j]; 01341 01342 } 01343 01344 } 01345 01346 // opserr<<"final "<<P<<endln; 01347 01348 return P; 01349 01350 } 01351 01352 01353 01354 int 01355 01356 NineFourNodeQuadUP::sendSelf(int commitTag, Channel &theChannel) 01357 01358 { 01359 01360 int res = 0; 01361 01362 01363 01364 // note: we don't check for dataTag == 0 for Element 01365 01366 // objects as that is taken care of in a commit by the Domain 01367 01368 // object - don't want to have to do the check if sending data 01369 01370 int dataTag = this->getDbTag(); 01371 01372 01373 01374 // Quad packs its data into a Vector and sends this to theChannel 01375 01376 // along with its dbTag and the commitTag passed in the arguments 01377 01378 static Vector data(10); 01379 01380 data(0) = this->getTag(); 01381 01382 data(1) = thickness; 01383 01384 data(2) = rho; 01385 01386 data(3) = b[0]; 01387 01388 data(4) = b[1]; 01389 01390 //data(5) = pressure; 01391 01392 data(6) = kc; 01393 01394 data(7) = perm[0]; 01395 01396 data(8) = perm[1]; 01397 01398 data(9) = 9; 01399 01400 01401 01402 01403 01404 res += theChannel.sendVector(dataTag, commitTag, data); 01405 01406 if (res < 0) { 01407 01408 opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send Vector\n"; 01409 01410 return res; 01411 01412 } 01413 01414 01415 01416 // Quad then sends the tags of its four end nodes 01417 01418 res += theChannel.sendID(dataTag, commitTag, connectedExternalNodes); 01419 01420 if (res < 0) { 01421 01422 opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send ID\n"; 01423 01424 return res; 01425 01426 } 01427 01428 01429 01430 // Now quad sends the ids of its materials 01431 01432 int matDbTag; 01433 01434 int numMats = 9; 01435 01436 ID classTags(2*numMats); 01437 01438 01439 01440 int i; 01441 01442 for (i = 0; i < 9; i++) { 01443 01444 classTags(i) = theMaterial[i]->getClassTag(); 01445 01446 matDbTag = theMaterial[i]->getDbTag(); 01447 01448 // NOTE: we do have to ensure that the material has a database 01449 01450 // tag if we are sending to a database channel. 01451 01452 if (matDbTag == 0) { 01453 01454 matDbTag = theChannel.getDbTag(); 01455 01456 if (matDbTag != 0) 01457 01458 theMaterial[i]->setDbTag(matDbTag); 01459 01460 } 01461 01462 classTags(i+numMats) = matDbTag; 01463 01464 } 01465 01466 01467 01468 res += theChannel.sendID(dataTag, commitTag, classTags); 01469 01470 if (res < 0) { 01471 01472 opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send ID\n"; 01473 01474 return res; 01475 01476 } 01477 01478 01479 01480 // Finally, quad asks its material objects to send themselves 01481 01482 for (i = 0; i < 9; i++) { 01483 01484 res += theMaterial[i]->sendSelf(commitTag, theChannel); 01485 01486 if (res < 0) { 01487 01488 opserr << "WARNING NineFourNodeQuadUP::sendSelf() - " << this->getTag() << " failed to send its Material\n"; 01489 01490 return res; 01491 01492 } 01493 01494 } 01495 01496 01497 01498 return res; 01499 01500 } 01501 01502 01503 01504 int 01505 01506 NineFourNodeQuadUP::recvSelf(int commitTag, Channel &theChannel, 01507 01508 FEM_ObjectBroker &theBroker) 01509 01510 { 01511 01512 int res = 0; 01513 01514 01515 01516 int dataTag = this->getDbTag(); 01517 01518 01519 01520 // Quad creates a Vector, receives the Vector and then sets the 01521 01522 // internal data with the data in the Vector 01523 01524 static Vector data(7); 01525 01526 res += theChannel.recvVector(dataTag, commitTag, data); 01527 01528 if (res < 0) { 01529 01530 opserr << "WARNING NineFourNodeQuadUP::recvSelf() - failed to receive Vector\n"; 01531 01532 return res; 01533 01534 } 01535 01536 01537 01538 this->setTag((int)data(0)); 01539 01540 thickness = data(1); 01541 01542 rho = data(2); 01543 01544 b[0] = data(3); 01545 01546 b[1] = data(4); 01547 01548 //pressure = data(5); 01549 01550 kc = data(6); 01551 01552 perm[0] = data(7); 01553 01554 perm[1] = data(8); 01555 01556 01557 01558 // Quad now receives the tags of its four external nodes 01559 01560 res += theChannel.recvID(dataTag, commitTag, connectedExternalNodes); 01561 01562 if (res < 0) { 01563 01564 opserr << "WARNING NineFourNodeQuadUP::recvSelf() - " << this->getTag() << " failed to receive ID\n"; 01565 01566 return res; 01567 01568 } 01569 01570 01571 01572 // Quad now receives the ids of its materials 01573 01574 int newOrder = (int)data(9); 01575 01576 int numMats = newOrder; 01577 01578 ID classTags(2*numMats); 01579 01580 01581 01582 res += theChannel.recvID(dataTag, commitTag, classTags); 01583 01584 if (res < 0) { 01585 01586 opserr << "NineFourNodeQuadUP::recvSelf() - failed to recv ID data\n"; 01587 01588 return res; 01589 01590 } 01591 01592 01593 01594 int i; 01595 01596 01597 01598 // If the number of materials (quadrature order) is not the same, 01599 01600 // delete the old materials, allocate new ones and then receive 01601 01602 if (9 != newOrder) { 01603 01604 // Delete the materials 01605 01606 for (i = 0; i < 9; i++) { 01607 01608 if (theMaterial[i]) 01609 01610 delete theMaterial[i]; 01611 01612 } 01613 01614 if (theMaterial) 01615 01616 delete [] theMaterial; 01617 01618 01619 01620 // Allocate new materials 01621 01622 theMaterial = new NDMaterial *[9]; 01623 01624 if (theMaterial == 0) { 01625 01626 opserr << "NineFourNodeQuadUP::recvSelf() - Could not allocate NDMaterial* array\n"; 01627 01628 return -1; 01629 01630 } 01631 01632 for (i = 0; i < 9; i++) { 01633 01634 int matClassTag = classTags(i); 01635 01636 int matDbTag = classTags(i+numMats); 01637 01638 // Allocate new material with the sent class tag 01639 01640 theMaterial[i] = theBroker.getNewNDMaterial(matClassTag); 01641 01642 if (theMaterial[i] == 0) { 01643 01644 opserr << "NineFourNodeQuadUP::recvSelf() - Broker could not create NDMaterial of class type" << matClassTag << endln; 01645 01646 return -1; 01647 01648 } 01649 01650 // Now receive materials into the newly allocated space 01651 01652 theMaterial[i]->setDbTag(matDbTag); 01653 01654 res += theMaterial[i]->recvSelf(commitTag, theChannel, theBroker); 01655 01656 if (res < 0) { 01657 01658 opserr << "NLBeamColumn3d::recvSelf() - material " << i << "failed to recv itself\n"; 01659 01660 return res; 01661 01662 } 01663 01664 } 01665 01666 } 01667 01668 // Number of materials is the same, receive materials into current space 01669 01670 else { 01671 01672 for (i = 0; i < 9; i++) { 01673 01674 int matClassTag = classTags(i); 01675 01676 int matDbTag = classTags(i+numMats); 01677 01678 // Check that material is of the right type; if not, 01679 01680 // delete it and create a new one of the right type 01681 01682 if (theMaterial[i]->getClassTag() != matClassTag) { 01683 01684 delete theMaterial[i]; 01685 01686 theMaterial[i] = theBroker.getNewNDMaterial(matClassTag); 01687 01688 if (theMaterial[i] == 0) { 01689 01690 opserr << "NineFourNodeQuadUP::recvSelf() - Broker could not create NDMaterial of class type " << matClassTag << endln; 01691 01692 exit(-1); 01693 01694 } 01695 01696 } 01697 01698 // Receive the material 01699 01700 theMaterial[i]->setDbTag(matDbTag); 01701 01702 res += theMaterial[i]->recvSelf(commitTag, theChannel, theBroker); 01703 01704 if (res < 0) { 01705 01706 opserr << "NineFourNodeQuadUP::recvSelf() - material " << i << "failed to recv itself\n"; 01707 01708 return res; 01709 01710 } 01711 01712 } 01713 01714 } 01715 01716 01717 01718 return res; 01719 01720 } 01721 01722 01723 01724 void 01725 01726 NineFourNodeQuadUP::Print(OPS_Stream &s, int flag) 01727 01728 { 01729 01730 s << "\nNineFourNodeQuadUP, element id: " << this->getTag() << endln; 01731 01732 s << "\tConnected external nodes: " << connectedExternalNodes; 01733 01734 s << "\tthickness: " << thickness << endln; 01735 01736 s << "\tmass density: " << rho << endln; 01737 01738 //s << "\tsurface pressure: " << pressure << endln; 01739 01740 s << "\tbody forces: " << b[0] << ' ' << b[1] << endln; 01741 01742 theMaterial[0]->Print(s,flag); 01743 01744 s << "\tStress (xx yy xy)" << endln; 01745 01746 for (int i = 0; i < 9; i++) 01747 01748 s << "\t\tGauss point " << i+1 << ": " << theMaterial[i]->getStress(); 01749 01750 } 01751 01752 01753 01754 int 01755 01756 NineFourNodeQuadUP::displaySelf(Renderer &theViewer, int displayMode, float fact) 01757 01758 { 01759 01760 // first set the quantity to be displayed at the nodes; 01761 01762 // if displayMode is 1 through 3 we will plot material stresses otherwise 0.0 01763 01764 01765 01766 static Vector values(9); 01767 01768 01769 01770 for (int j=0; j<9; j++) 01771 01772 values(j) = 0.0; 01773 01774 01775 01776 if (displayMode < 4 && displayMode > 0) { 01777 01778 for (int i=0; i<9; i++) { 01779 01780 const Vector &stress = theMaterial[i]->getStress(); 01781 01782 values(i) = stress(displayMode-1); 01783 01784 } 01785 01786 } 01787 01788 01789 01790 // now determine the end points of the quad based on 01791 01792 // the display factor (a measure of the distorted image) 01793 01794 // store this information in 4 3d vectors v1 through v4 01795 01796 /*const Vector &end1Crd = nd1Ptr->getCrds(); 01797 01798 const Vector &end2Crd = nd2Ptr->getCrds(); 01799 01800 const Vector &end3Crd = nd3Ptr->getCrds(); 01801 01802 const Vector &end4Crd = nd4Ptr->getCrds(); 01803 01804 01805 01806 const Vector &end1Disp = nd1Ptr->getDisp(); 01807 01808 const Vector &end2Disp = nd2Ptr->getDisp(); 01809 01810 const Vector &end3Disp = nd3Ptr->getDisp(); 01811 01812 const Vector &end4Disp = nd4Ptr->getDisp(); 01813 01814 01815 01816 static Matrix coords(4,3); 01817 01818 01819 01820 for (int i = 0; i < 2; i++) { 01821 01822 coords(0,i) = end1Crd(i) + end1Disp(i)*fact; 01823 01824 coords(1,i) = end2Crd(i) + end2Disp(i)*fact; 01825 01826 coords(2,i) = end3Crd(i) + end3Disp(i)*fact; 01827 01828 coords(3,i) = end4Crd(i) + end4Disp(i)*fact; 01829 01830 } 01831 01832 01833 01834 int error = 0; 01835 01836 01837 01838 // finally we draw the element using drawPolygon 01839 01840 error += theViewer.drawPolygon (coords, values); 01841 01842 01843 01844 return error;*/ 01845 01846 return 0; 01847 01848 } 01849 01850 01851 01852 Response* 01853 NineFourNodeQuadUP::setResponse(const char **argv, int argc, Information &eleInfo, OPS_Stream &output) 01854 01855 { 01856 Response *theResponse = 0; 01857 01858 char outputData[32]; 01859 01860 output.tag("ElementOutput"); 01861 output.attr("eleType","NineFOurNodeQuadUP"); 01862 output.attr("eleTag",this->getTag()); 01863 for (int i=1; i<=9; i++) { 01864 sprintf(outputData,"node%d",i); 01865 output.attr(outputData, theNodes[i-1]->getTag()); 01866 } 01867 01868 if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0) { 01869 01870 01871 for (int i=1; i<=9; i++) { 01872 sprintf(outputData,"P1_%d",i); 01873 output.tag("ResponseType",outputData); 01874 sprintf(outputData,"P2_%d",i); 01875 output.tag("ResponseType",outputData); 01876 if (i <= nenp) { 01877 sprintf(outputData,"Pp_%d",i); 01878 output.tag("ResponseType",outputData); 01879 } 01880 } 01881 01882 theResponse = new ElementResponse(this, 1, P); 01883 01884 } else if (strcmp(argv[0],"stiff") == 0 || strcmp(argv[0],"stiffness") == 0) { 01885 theResponse = new ElementResponse(this, 2, K); 01886 01887 01888 } else if (strcmp(argv[0],"mass") == 0) { 01889 01890 theResponse = new ElementResponse(this, 3, K); 01891 01892 01893 01894 } else if (strcmp(argv[0],"damp") == 0) { 01895 01896 theResponse = new ElementResponse(this, 4, K); 01897 01898 01899 01900 } else if (strcmp(argv[0],"material") == 0 || strcmp(argv[0],"integrPoint") == 0) { 01901 01902 int pointNum = atoi(argv[1]); 01903 01904 if (pointNum > 0 && pointNum <= nenu) { 01905 01906 output.tag("GaussPoint"); 01907 output.attr("number",pointNum); 01908 01909 theResponse = theMaterial[pointNum-1]->setResponse(&argv[2], argc-2, eleInfo, output); 01910 01911 output.endTag(); // GaussPoint 01912 } 01913 } 01914 01915 output.endTag(); // ElementOutput 01916 return theResponse; 01917 } 01918 01919 01920 01921 int 01922 01923 NineFourNodeQuadUP::getResponse(int responseID, Information &eleInfo) 01924 01925 { 01926 01927 switch (responseID) { 01928 01929 01930 01931 case 1: 01932 01933 return eleInfo.setVector(this->getResistingForce()); 01934 01935 01936 01937 case 2: 01938 01939 return eleInfo.setMatrix(this->getTangentStiff()); 01940 01941 01942 01943 case 3: 01944 01945 return eleInfo.setMatrix(this->getMass()); 01946 01947 01948 01949 case 4: 01950 01951 return eleInfo.setMatrix(this->getDamp()); 01952 01953 01954 01955 default: 01956 01957 return -1; 01958 01959 } 01960 01961 } 01962 01963 int 01964 NineFourNodeQuadUP::setParameter(const char **argv, int argc, Parameter ¶m) 01965 { 01966 if (argc < 1) 01967 return -1; 01968 01969 // quad mass density per unit volume 01970 if (strcmp(argv[0],"rho") == 0) 01971 return param.addObject(1, this); 01972 01973 // a material parameter 01974 else if (strstr(argv[0],"material") != 0) { 01975 01976 if (argc < 3) 01977 return -1; 01978 01979 int pointNum = atoi(argv[1]); 01980 if (pointNum > 0 && pointNum <= nenu) 01981 return theMaterial[pointNum-1]->setParameter(&argv[2], argc-2, param); 01982 else 01983 return -1; 01984 } 01985 01986 // otherwise parameter is unknown for the Truss class 01987 else 01988 return -1; 01989 } 01990 01991 int 01992 NineFourNodeQuadUP::updateParameter(int parameterID, Information &info) 01993 { 01994 switch (parameterID) { 01995 case 1: 01996 rho = info.theDouble; 01997 this->getMass(); // update mass matrix 01998 return 0; 01999 02000 default: 02001 return -1; 02002 } 02003 } 02004 02005 void NineFourNodeQuadUP::globalShapeFunction(double *dvol, double *w, int nint, int nen, int mode) 02006 02007 { 02008 02009 static double coord[2][9], xs[2][2], det, temp; 02010 02011 int i, j, k, m; 02012 02013 02014 02015 for (i=0; i<3; i++) { 02016 02017 for (j=0; j<nen; j++) { 02018 02019 for (k=0; k<nint; k++) { 02020 02021 if (mode==0) shgu[i][j][k] = shlu[i][j][k]; 02022 02023 if (mode==1) shgp[i][j][k] = shlp[i][j][k]; 02024 02025 if (mode==2) shgq[i][j][k] = shlq[i][j][k]; 02026 02027 } 02028 02029 } 02030 02031 } 02032 02033 02034 02035 for (i=0; i<nen; i++) { 02036 02037 const Vector &coordRef = theNodes[i]->getCrds(); 02038 02039 coord[0][i] = coordRef(0); 02040 02041 coord[1][i] = coordRef(1); 02042 02043 } 02044 02045 02046 02047 for (m=0; m<nint; m++) { 02048 02049 for (i=0; i<2; i++) { 02050 02051 for (j=0; j<2; j++) { 02052 02053 xs[i][j] = 0.0; 02054 02055 for (k=0; k<nen; k++) { 02056 02057 if (mode==0) xs[i][j] += coord[j][k]*shgu[i][k][m]; 02058 02059 if (mode==1) xs[i][j] += coord[j][k]*shgp[i][k][m]; 02060 02061 if (mode==2) xs[i][j] += coord[j][k]*shgq[i][k][m]; 02062 02063 } 02064 02065 } 02066 02067 } 02068 02069 02070 02071 det = xs[0][0]*xs[1][1] - xs[0][1]*xs[1][0]; 02072 02073 02074 02075 if (det < 0.0) { 02076 02077 opserr << "WARNING NineFourNodeQuadUP: Determinant<=0 in tag " 02078 02079 << this->getTag(); 02080 02081 exit(-1); 02082 02083 } 02084 02085 02086 02087 for (i=0; i<nen; i++) { 02088 02089 if (mode==0) { 02090 02091 temp = (shgu[0][i][m]*xs[1][1] - shgu[1][i][m]*xs[0][1])/det; 02092 02093 shgu[1][i][m] = (-shgu[0][i][m]*xs[1][0] + shgu[1][i][m]*xs[0][0])/det; 02094 02095 shgu[0][i][m] = temp; 02096 02097 } 02098 02099 if (mode==1) { 02100 02101 temp = (shgp[0][i][m]*xs[1][1] - shgp[1][i][m]*xs[0][1])/det; 02102 02103 shgp[1][i][m] = (-shgp[0][i][m]*xs[1][0] + shgp[1][i][m]*xs[0][0])/det; 02104 02105 shgp[0][i][m] = temp; 02106 02107 } 02108 02109 if (mode==2) { 02110 02111 temp = (shgq[0][i][m]*xs[1][1] - shgq[1][i][m]*xs[0][1])/det; 02112 02113 shgq[1][i][m] = (-shgq[0][i][m]*xs[1][0] + shgq[1][i][m]*xs[0][0])/det; 02114 02115 shgq[0][i][m] = temp; 02116 02117 } 02118 02119 } 02120 02121 02122 02123 dvol[m] = w[m]*thickness*det; 02124 02125 02126 02127 } //end of m loop 02128 02129 02130 02131 } 02132 02133 02134 02135 02136 02137 void NineFourNodeQuadUP::shapeFunction(double *w, int nint, int nen, int mode) 02138 02139 { 02140 02141 static const double ra[] = {-0.5,0.5,0.5,-0.5,0.,0.5,0.,-0.5,0.}; 02142 02143 static const double sa[] = {-0.5,-0.5,0.5,0.5,-0.5,0.,0.5,0.,0.}; 02144 02145 02146 02147 double g, r, s, shl19, shl29, shl39, tempr, temps; 02148 02149 int ic, il, is; 02150 02151 02152 02153 g = 0.; 02154 02155 if (nint == 4) { 02156 02157 g=2./sqrt(3.0); 02158 02159 w[0] = w[1] = w[2] = w[3] = 1.; 02160 02161 } 02162 02163 if (nint == 9) { 02164 02165 g=2.*sqrt(3.0/5.0); 02166 02167 w[0] = w[1] = w[2] = w[3] = 25./81.; 02168 02169 w[4] = w[5] = w[6] = w[7] = 40./81.; 02170 02171 w[8] = 64./81.; 02172 02173 } 02174 02175 02176 02177 for (int i=0; i<nint; i++) { 02178 02179 r = g*ra[i]; 02180 02181 s = g*sa[i]; 02182 02183 shl19 = shl29 = shl39 = 0.; 02184 02185 if (nen > 4) { 02186 02187 tempr = 1.-r*r; 02188 02189 temps = 1.-s*s; 02190 02191 if (nen == 9) { 02192 02193 if (mode==0) { 02194 02195 shlu[0][8][i] = -2.*r*temps; 02196 02197 shl19 = 0.5*shlu[0][8][i]; 02198 02199 shlu[1][8][i] = -2.*s*tempr; 02200 02201 shl29 = 0.5*shlu[1][8][i]; 02202 02203 shlu[2][8][i] = temps*tempr; 02204 02205 shl39 = 0.5*shlu[2][8][i]; 02206 02207 } 02208 02209 if (mode==2) { 02210 02211 shlq[0][8][i] = -2.*r*temps; 02212 02213 shl19 = 0.5*shlq[0][8][i]; 02214 02215 shlq[1][8][i] = -2.*s*tempr; 02216 02217 shl29 = 0.5*shlq[1][8][i]; 02218 02219 shlq[2][8][i] = temps*tempr; 02220 02221 shl39 = 0.5*shlq[2][8][i]; 02222 02223 } 02224 02225 } 02226 02227 if (mode==0) { 02228 02229 shlu[0][4][i] = -r*(1.-s) - shl19; 02230 02231 shlu[1][4][i] = -0.5*tempr - shl29; 02232 02233 shlu[2][4][i] = 0.5*tempr*(1.-s) - shl39; 02234 02235 shlu[0][5][i] = 0.5*temps - shl19; 02236 02237 shlu[1][5][i] = -s*(1.+r) - shl29; 02238 02239 shlu[2][5][i] = 0.5*temps*(1.+r) - shl39; 02240 02241 shlu[0][6][i] = -r*(1.+s) - shl19; 02242 02243 shlu[1][6][i] = 0.5*tempr - shl29; 02244 02245 shlu[2][6][i] = 0.5*tempr*(1.+s) - shl39; 02246 02247 shlu[0][7][i] = -0.5*temps - shl19; 02248 02249 shlu[1][7][i] = -s*(1.-r) - shl29; 02250 02251 shlu[2][7][i] = 0.5*temps*(1.-r) - shl39; 02252 02253 } 02254 02255 if (mode==2) { 02256 02257 shlq[0][4][i] = -r*(1.-s) - shl19; 02258 02259 shlq[1][4][i] = -0.5*tempr - shl29; 02260 02261 shlq[2][4][i] = 0.5*tempr*(1.-s) - shl39; 02262 02263 shlq[0][5][i] = 0.5*temps - shl19; 02264 02265 shlq[1][5][i] = -s*(1.+r) - shl29; 02266 02267 shlq[2][5][i] = 0.5*temps*(1.+r) - shl39; 02268 02269 shlq[0][6][i] = -r*(1.+s) - shl19; 02270 02271 shlq[1][6][i] = 0.5*tempr - shl29; 02272 02273 shlq[2][6][i] = 0.5*tempr*(1.+s) - shl39; 02274 02275 shlq[0][7][i] = -0.5*temps - shl19; 02276 02277 shlq[1][7][i] = -s*(1.-r) - shl29; 02278 02279 shlq[2][7][i] = 0.5*temps*(1.-r) - shl39; 02280 02281 } 02282 02283 } 02284 02285 02286 02287 for (int k=0; k<4; k++) { 02288 02289 tempr = 0.5 + ra[k]*r; 02290 02291 temps = 0.5 + sa[k]*s; 02292 02293 if (mode==0) { 02294 02295 shlu[0][k][i] = ra[k]*temps - 0.5*shl19; 02296 02297 shlu[1][k][i] = tempr*sa[k] - 0.5*shl29; 02298 02299 shlu[2][k][i] = tempr*temps - 0.5*shl39; 02300 02301 } 02302 02303 if (mode==1) { 02304 02305 shlp[0][k][i] = ra[k]*temps - 0.5*shl19; 02306 02307 shlp[1][k][i] = tempr*sa[k] - 0.5*shl29; 02308 02309 shlp[2][k][i] = tempr*temps - 0.5*shl39; 02310 02311 } 02312 02313 if (mode==2) { 02314 02315 shlq[0][k][i] = ra[k]*temps - 0.5*shl19; 02316 02317 shlq[1][k][i] = tempr*sa[k] - 0.5*shl29; 02318 02319 shlq[2][k][i] = tempr*temps - 0.5*shl39; 02320 02321 } 02322 02323 } 02324 02325 02326 02327 if (nen > 4) { 02328 02329 for (int m=4; m<8; m++) { 02330 02331 ic = m - 4; 02332 02333 il = m - 3; 02334 02335 is = 1; 02336 02337 if (m==7) { 02338 02339 ic = 0; 02340 02341 il = 3; 02342 02343 is = 3; 02344 02345 } 02346 02347 for (int j=ic; j<=il; j+=is) { 02348 02349 if (mode==0) { 02350 02351 shlu[0][j][i] = shlu[0][j][i] - 0.5*shlu[0][m][i]; 02352 02353 shlu[1][j][i] = shlu[1][j][i] - 0.5*shlu[1][m][i]; 02354 02355 shlu[2][j][i] = shlu[2][j][i] - 0.5*shlu[2][m][i]; 02356 02357 } 02358 02359 if (mode==2) { 02360 02361 shlq[0][j][i] = shlq[0][j][i] - 0.5*shlq[0][m][i]; 02362 02363 shlq[1][j][i] = shlq[1][j][i] - 0.5*shlq[1][m][i]; 02364 02365 shlq[2][j][i] = shlq[2][j][i] - 0.5*shlq[2][m][i]; 02366 02367 } 02368 02369 } 02370 02371 } //end of m loop 02372 02373 } 02374 02375 } //end of i loop 02376 02377 } 02378 02379 02380 02381 02382 02383 double NineFourNodeQuadUP::mixtureRho(int i) 02384 02385 { 02386 02387 double rhoi, e, n; 02388 02389 rhoi= theMaterial[i]->getRho(); 02390 02391 e = 0.7; //theMaterial[i]->getVoidRatio(); 02392 02393 n = e / (1.0 + e); 02394 02395 //return n * rho + (1.0-n) * rhoi; 02396 02397 return rhoi; 02398 02399 }
Generated on Mon Oct 23 15:05:10 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents