EightNodeBrick_u_p_U.cppGo to the documentation of this file.00001 00002 // 00003 // COPYRIGHT (C): :-)) 00004 // PROJECT: Object Oriented Finite Element Program 00005 // FILE: EightNodeBrick_u_p_U.cpp 00006 // CLASS: EightNodeBrick_u_p_U 00007 // MEMBER FUNCTIONS: 00008 // 00009 // MEMBER VARIABLES 00010 // 00011 // PURPOSE: Finite Element Class for coupled system 00012 // "Coupled system": Solid and fluid coexist. 00013 // u-- Solid displacement 00014 // p-- Pore pressure 00015 // U-- Absolute fluid displacement 00016 // 00017 // RETURN: 00018 // VERSION: 00019 // LANGUAGE: C++ 00020 // TARGET OS: DOS || UNIX || . . . 00021 // DESIGNER: Boris Jeremic, Zhao Cheng 00022 // PROGRAMMER: Boris Jeremic, Zhaohui Yang, Xiaoyan Wu, Zhao Cheng 00023 // DATE: Aug. 2001 00024 // UPDATE HISTORY: Modified from EightNodeBrick.cpp reorganized a lot by Xiaoyan 00025 // 01/24/2002 Xiaoyan 00026 // Add the permeability tensor and ks, kf to the constructor Xiaoyan 00027 // 00028 // Clean-up and re-write by Zhao Cheng, 10/20/2004 00029 // 00030 // Fixed a bug, and some small modification, ZC 05/2006 00032 // 00033 #ifndef EIGHTNODEBRICK_U_P_U_CPP 00034 #define EIGHTNODEBRICK_U_P_U_CPP 00035 00036 #include <EightNodeBrick_u_p_U.h> 00037 00038 const int EightNodeBrick_u_p_U::Num_IntegrationPts = 2; 00039 const int EightNodeBrick_u_p_U::Num_TotalGaussPts = 8; 00040 const int EightNodeBrick_u_p_U::Num_Nodes = 8; 00041 const int EightNodeBrick_u_p_U::Num_Dim = 3; 00042 const int EightNodeBrick_u_p_U::Num_Dof = 7; 00043 const int EightNodeBrick_u_p_U::Num_ElemDof = 56; 00044 const double EightNodeBrick_u_p_U::pts[2] = {-0.577350269189626, +0.577350269189626}; 00045 const double EightNodeBrick_u_p_U::wts[2] = {1.0, 1.0}; 00046 Matrix EightNodeBrick_u_p_U::MCK(Num_ElemDof, Num_ElemDof); 00047 Vector EightNodeBrick_u_p_U::P(Num_ElemDof); 00048 00049 //====================================================================== 00050 EightNodeBrick_u_p_U::EightNodeBrick_u_p_U(int element_number, 00051 int node_numb_1, 00052 int node_numb_2, 00053 int node_numb_3, 00054 int node_numb_4, 00055 int node_numb_5, 00056 int node_numb_6, 00057 int node_numb_7, 00058 int node_numb_8, 00059 NDMaterial *Globalmmodel, 00060 double b1, 00061 double b2, 00062 double b3, 00063 double nn, 00064 double alf, 00065 double rs, 00066 double rf, 00067 double permb_x, 00068 double permb_y, 00069 double permb_z, 00070 double kks, 00071 double kkf, 00072 double pp) 00073 : Element(element_number, 00074 ELE_TAG_EightNodeBrick_u_p_U ), 00075 connectedExternalNodes(Num_Nodes), 00076 perm(Num_Dim), 00077 bf(Num_Dim), 00078 poro(nn), 00079 alpha(alf), 00080 rho_s(rs), 00081 rho_f(rf), 00082 ks(kks), 00083 kf(kkf), 00084 pressure(pp), 00085 Q(0), 00086 Ki(0) 00087 { 00088 // permeability 00089 perm(0) = permb_x; 00090 perm(1) = permb_y; 00091 perm(2) = permb_z; 00092 00093 if (perm(0)==0.0 || perm(1)==0.0 || perm(2)==0.0) { 00094 opserr<<" Error, EightNodeBrick_u_p_U:: permeability (kx/ky/kz) is zero! \n"; 00095 exit(-1); 00096 } 00097 00098 // body forces 00099 bf(0) = b1; 00100 bf(1) = b2; 00101 bf(2) = b3; 00102 00103 connectedExternalNodes( 0) = node_numb_1; 00104 connectedExternalNodes( 1) = node_numb_2; 00105 connectedExternalNodes( 2) = node_numb_3; 00106 connectedExternalNodes( 3) = node_numb_4; 00107 connectedExternalNodes( 4) = node_numb_5; 00108 connectedExternalNodes( 5) = node_numb_6; 00109 connectedExternalNodes( 6) = node_numb_7; 00110 connectedExternalNodes( 7) = node_numb_8; 00111 00112 theMaterial = new NDMaterial *[Num_TotalGaussPts]; 00113 if (theMaterial == 0) { 00114 opserr<<" EightNodeBrick_u_p_U::EightNodeBrick_u_p_U -- failed allocate material model pointer\n"; 00115 exit(-1); 00116 } 00117 for (int i=0; i<Num_TotalGaussPts; i++) { 00118 theMaterial[i] = Globalmmodel->getCopy(); 00119 if (theMaterial[i] == 0) { 00120 opserr<<" EightNodeBrick_u_p_U::EightNodeBrick_u_p_U -- failed allocate material model pointer\n"; 00121 exit(-1); 00122 } 00123 } 00124 00125 } 00126 00127 //====================================================================== 00128 EightNodeBrick_u_p_U::EightNodeBrick_u_p_U () 00129 : Element(0, ELE_TAG_EightNodeBrick_u_p_U ), 00130 connectedExternalNodes(Num_Nodes), perm(Num_Dim), bf(Num_Dim), 00131 poro(0.0), alpha(1.0), rho_s(0.0),rho_f(0.0), ks(0.0), kf(0.0), pressure(0.0), Q(0), Ki(0) 00132 { 00133 theMaterial = 0; 00134 00135 for (int j=0; j<Num_Nodes; j++) 00136 theNodes[j] = 0; 00137 } 00138 00139 //====================================================================== 00140 EightNodeBrick_u_p_U::~EightNodeBrick_u_p_U () 00141 { 00142 if (theMaterial) 00143 delete [] theMaterial; 00144 00145 for (int j=0; j<Num_Nodes; j++) 00146 theNodes[j] = 0; 00147 00148 if (Q != 0) 00149 delete Q; 00150 00151 if (Ki != 0) 00152 delete Ki; 00153 00154 } 00155 00156 //====================================================================== 00157 int EightNodeBrick_u_p_U::getNumExternalNodes (void) const 00158 { 00159 return Num_Nodes; 00160 } 00161 00162 //====================================================================== 00163 const ID& EightNodeBrick_u_p_U::getExternalNodes (void) 00164 { 00165 return connectedExternalNodes; 00166 } 00167 00168 //====================================================================== 00169 Node ** EightNodeBrick_u_p_U::getNodePtrs (void) 00170 { 00171 return theNodes; 00172 } 00173 00174 //====================================================================== 00175 int EightNodeBrick_u_p_U::getNumDOF (void) 00176 { 00177 return Num_ElemDof; 00178 } 00179 00180 //====================================================================== 00181 void EightNodeBrick_u_p_U::setDomain (Domain *theDomain) 00182 { 00183 int i, Ndof; 00184 00185 if (theDomain == 0) { 00186 for (i=0; i<Num_Nodes; i++) { 00187 theNodes[i] = 0; 00188 } 00189 } 00190 00191 for (i=0; i<Num_Nodes; i++) { 00192 theNodes[i] = theDomain->getNode(connectedExternalNodes(i)); 00193 if (theNodes[i] == 0) { 00194 opserr << "Error EightNodeBrick_u_p_U : node not found in the domain" << "\n"; 00195 return ; 00196 } 00197 Ndof = theNodes[i]->getNumberDOF(); 00198 if( Ndof != Num_Dof) { 00199 opserr << "Error EightNodeBrick_u_p_U : has wrong number of DOFs at its nodes" << "\n"; 00200 return ; 00201 } 00202 } 00203 00204 this->DomainComponent::setDomain(theDomain); 00205 00206 } 00207 00208 //====================================================================== 00209 int EightNodeBrick_u_p_U::commitState (void) 00210 { 00211 int retVal = 0; 00212 int i; 00213 00214 if ((retVal = this->Element::commitState()) != 0) { 00215 opserr << "EightNodeBrick-u_p_U::commitState () - failed in base class"; 00216 return (-1); 00217 } 00218 00219 for (i=0; i<Num_TotalGaussPts; i++ ) 00220 retVal += theMaterial[i]->commitState(); 00221 00222 return retVal; 00223 } 00224 00225 //====================================================================== 00226 int EightNodeBrick_u_p_U::revertToLastCommit (void) 00227 { 00228 int retVal = 0; 00229 int i; 00230 00231 for (i=0; i<Num_TotalGaussPts; i++ ) 00232 retVal += theMaterial[i]->revertToLastCommit() ; 00233 00234 return retVal; 00235 } 00236 00237 //====================================================================== 00238 int EightNodeBrick_u_p_U::revertToStart (void) 00239 { 00240 int retVal = 0; 00241 int i; 00242 00243 for (i=0; i<Num_TotalGaussPts; i++ ) 00244 retVal += theMaterial[i]->revertToStart() ; 00245 00246 return retVal; 00247 } 00248 00249 //====================================================================== 00250 const Matrix& EightNodeBrick_u_p_U::getTangentStiff (void) 00251 { 00252 return this->getStiff(1); 00253 } 00254 00255 //====================================================================== 00256 const Matrix& EightNodeBrick_u_p_U::getInitialStiff (void) 00257 { 00258 return this->getStiff(0); 00259 } 00260 00261 //====================================================================== 00262 const Matrix& EightNodeBrick_u_p_U::getDamp (void) 00263 { 00264 MCK.Zero(); // necessary 00265 00266 tensor tC = getDampTensorC123(); 00267 00268 int i, j, m, n; 00269 00270 double Ctemp = 0.0; 00271 tensor CRm; 00272 tensor CRk; 00273 00274 if (alphaM != 0.0) 00275 CRm = getMassTensorMsf(); 00276 00277 if (betaK != 0.0) 00278 CRk = getStiffnessTensorKep(); 00279 00280 if (betaK0 != 0.0 || betaKc != 0.0) { 00281 opserr << "Warning: EightNodeBrick-u_p_U:: betaK0 or betaKc are not used" << "\n"; 00282 } 00283 00284 for ( i=0 ; i<Num_Nodes; i++ ) { 00285 for ( j=0; j<Num_Nodes; j++ ) { 00286 for( m=0; m<Num_Dim; m++) { 00287 for( n=0; n<Num_Dim; n++) 00288 { 00289 Ctemp = tC.cval(i+1, m+1, n+1, j+1) *(poro*poro); 00290 00291 //C1 00292 MCK(i*Num_Dof+m, j*Num_Dof+n) = Ctemp; 00293 00294 if (alphaM != 0.0) 00295 MCK(i*Num_Dof+m, j*Num_Dof+n) += CRm.cval(i+1, j+1) *((1.0-poro)*rho_s *alphaM); 00296 00297 if (betaK != 0.0) 00298 MCK(i*Num_Dof+m, j*Num_Dof+n) += CRk.cval(i+1, m+1, n+1, j+1) * betaK; 00299 00300 //C3 00301 MCK(i*Num_Dof+m+4, j*Num_Dof+n+4) = Ctemp; 00302 00303 if (alphaM != 0.0) 00304 MCK(i*Num_Dof+m+4, j*Num_Dof+n+4) += CRm.cval(i+1, j+1) *(poro*rho_f *alphaM); 00305 00306 //C2 and C2^T 00307 MCK(i*Num_Dof+m, j*Num_Dof+n+4) = - Ctemp; 00308 MCK(j*Num_Dof+n+4, i*Num_Dof+m) = - Ctemp; 00309 } 00310 } 00311 } 00312 } 00313 00314 return MCK; 00315 } 00316 00317 //====================================================================== 00318 const Matrix& EightNodeBrick_u_p_U::getMass (void) 00319 { 00320 MCK.Zero(); // necessary 00321 00322 tensor tM = getMassTensorMsf(); 00323 00324 int i, j; 00325 double Mtemp1 = 0.0; 00326 double Mtemp2 = 0.0; 00327 00328 for ( i=0 ; i<Num_Nodes; i++ ) { 00329 for ( j=0; j<Num_Nodes; j++ ) { 00330 00331 //Ms, Note *(1.0-poro)*rho_s here! 00332 Mtemp1 = tM.cval(i+1, j+1) *(1.0-poro)*rho_s; 00333 00334 MCK(i*Num_Dof+0, j*Num_Dof+0) = Mtemp1; 00335 MCK(i*Num_Dof+1, j*Num_Dof+1) = Mtemp1; 00336 MCK(i*Num_Dof+2, j*Num_Dof+2) = Mtemp1; 00337 00338 //Mf, Note *poro*rho_f here! 00339 Mtemp2 = tM.cval(i+1, j+1) *poro*rho_f; 00340 00341 MCK(i*Num_Dof+4, j*Num_Dof+4) = Mtemp2; 00342 MCK(i*Num_Dof+5, j*Num_Dof+5) = Mtemp2; 00343 MCK(i*Num_Dof+6, j*Num_Dof+6) = Mtemp2; 00344 } 00345 } 00346 00347 return MCK; 00348 } 00349 00350 //====================================================================== 00351 void EightNodeBrick_u_p_U::zeroLoad() 00352 { 00353 if ( Q != 0 ) 00354 Q->Zero(); 00355 } 00356 00357 //====================================================================== 00358 int EightNodeBrick_u_p_U::addLoad(ElementalLoad *theLoad, double loadFactor) 00359 { 00360 int type; 00361 const Vector& data = theLoad->getData(type, loadFactor); 00362 00363 if ( type == LOAD_TAG_BrickSelfWeight ) { 00364 Vector Fbody = this->getBodyForce(); 00365 if ( Q == 0 ) 00366 Q = new Vector(Num_ElemDof); 00367 *Q = Fbody * loadFactor; 00368 } 00369 00370 else { 00371 opserr << "EightNodeBrick_u_p_U::addLoad() " << this->getTag() << ", load type unknown\n"; 00372 return -1; 00373 } 00374 00375 return 0; 00376 } 00377 00378 //====================================================================== 00379 int EightNodeBrick_u_p_U::addInertiaLoadToUnbalance(const Vector &accel) 00380 { 00381 static Vector avu(Num_ElemDof); 00382 00383 int i, ik; 00384 00385 for (i=0; i<Num_Nodes; i++) { 00386 const Vector &RA = theNodes[i]->getRV(accel); 00387 00388 if ( RA.Size() != Num_Dof ) { 00389 opserr << "EightNodeBrick_u_p_U::addInertiaLoadToUnbalance matrix and vector sizes are incompatable\n"; 00390 return (-1); 00391 } 00392 00393 ik = i*Num_Dof; 00394 00395 avu(ik +0) = RA(0); 00396 avu(ik +1) = RA(1); 00397 avu(ik +2) = RA(2); 00398 avu(ik +3) = 0.0; 00399 avu(ik +4) = RA(4); 00400 avu(ik +5) = RA(5); 00401 avu(ik +6) = RA(6); 00402 } 00403 00404 if (Q == 0) 00405 Q = new Vector(Num_ElemDof); 00406 00407 this->getMass(); 00408 Q->addMatrixVector(1.0, MCK, avu, -1.0); 00409 00410 return 0; 00411 } 00412 00413 //======================================================================== 00414 const Vector& EightNodeBrick_u_p_U::getResistingForce () 00415 { 00416 static Vector avu(Num_ElemDof); 00417 00418 P.Zero(); 00419 00420 // Internal 00421 P = this->getInternalForce(); 00422 00423 // + K0*u 00424 int i, j; 00425 for (i=0; i<Num_Nodes; i++) { 00426 const Vector &disp = theNodes[i]->getTrialDisp(); 00427 if ( disp.Size() != Num_Dof ) { 00428 opserr << "EightNode_Brick_u_p_U::getResistingForce(): matrix and vector sizes are incompatable \n"; 00429 exit(-1); 00430 } 00431 for (j=0; j<Num_Dof; j++) { 00432 avu(i*Num_Dof +j) = disp(j); 00433 } 00434 } 00435 00436 this->getStiff00(); 00437 P.addMatrixVector(1.0, MCK, avu, 1.0); 00438 00439 if (Q != 0) 00440 P.addVector(1.0, *Q, -1.0); 00441 00442 return P; 00443 } 00444 00445 //======================================================================== 00446 const Vector& EightNodeBrick_u_p_U::getResistingForceIncInertia () 00447 { 00448 static Vector avu(Num_ElemDof); 00449 00450 this->getResistingForce(); 00451 00452 // + M*a 00453 int i, j; 00454 for (i=0; i<Num_Nodes; i++) { 00455 const Vector &acc = theNodes[i]->getTrialAccel(); 00456 if ( acc.Size() != Num_Dof ) { 00457 opserr << "EightNode_Brick_u_p_U::getResistingForceIncInertia matrix and vector sizes are incompatable \n"; 00458 exit(-1); 00459 } 00460 for (j=0; j<Num_Dof; j++) { 00461 avu(i*Num_Dof +j) = acc(j); 00462 } 00463 } 00464 00465 this->getMass(); 00466 P.addMatrixVector(1.0, MCK, avu, 1.0); 00467 00468 // + C*v 00469 for (i=0; i<Num_Nodes; i++) { 00470 const Vector &vel = theNodes[i]->getTrialVel(); 00471 if ( vel.Size() != Num_Dof ) { 00472 opserr << "EightNode_Brick_u_p_U::getResistingForceIncInertia matrix and vector sizes are incompatable \n"; 00473 exit(-1); 00474 } 00475 for (j=0; j<Num_Dof; j++) { 00476 avu(i*Num_Dof +j) = vel(j); 00477 } 00478 } 00479 00480 this->getDamp(); 00481 P.addMatrixVector(1.0, MCK, avu, 1.0); 00482 00483 return P; 00484 } 00485 00486 //============================================================================= 00487 int EightNodeBrick_u_p_U::sendSelf (int commitTag, Channel &theChannel) 00488 { 00489 // Not implemented yet 00490 return 0; 00491 } 00492 00493 //============================================================================= 00494 int EightNodeBrick_u_p_U::recvSelf (int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) 00495 { 00496 // Not implemented yet 00497 return 0; 00498 } 00499 00500 //============================================================================= 00501 int EightNodeBrick_u_p_U::displaySelf (Renderer &theViewer, int displayMode, float fact) 00502 { 00503 // Not implemented yet 00504 return 0; 00505 } 00506 00507 00508 //============================================================================= 00509 Response* EightNodeBrick_u_p_U::setResponse(const char **argv, int argc, Information &eleInfo, OPS_Stream &output) 00510 { 00511 Response *theResponse = 0; 00512 00513 char outputData[32]; 00514 00515 output.tag("ElementOutput"); 00516 output.attr("eleType","EightNodeBrick_u_p_U"); 00517 output.attr("eleTag",this->getTag()); 00518 for (int i=1; i<=Num_Nodes; i++) { 00519 sprintf(outputData,"node%d",i); 00520 output.attr(outputData,connectedExternalNodes[i-1]); 00521 } 00522 00523 if (strcmp(argv[0],"force") == 0 || strcmp(argv[0],"forces") == 0) { 00524 for (int i=1; i<=Num_Nodes; i++) 00525 for (int j=1; j<=Num_Dof; j++) { 00526 sprintf(outputData,"P%d_%d",j,i); 00527 output.tag("ResponseType",outputData); 00528 } 00529 00530 theResponse = new ElementResponse(this, 1, this->getResistingForce()); 00531 00532 } else if (strcmp(argv[0],"material") == 0 || strcmp(argv[0],"integrPoint") == 0) { 00533 int pointNum = atoi(argv[1]); 00534 if (pointNum > 0 && pointNum <= Num_TotalGaussPts) { 00535 00536 output.tag("GaussPoint"); 00537 output.attr("number",pointNum); 00538 00539 theResponse = theMaterial[pointNum-1]->setResponse(&argv[2], argc-2, eleInfo, output); 00540 00541 output.endTag(); // GaussPoint 00542 } 00543 00544 } else if (strcmp(argv[0],"stresses") ==0) { 00545 00546 for (int i=0; i<Num_TotalGaussPts; i++) { 00547 output.tag("GaussPoint"); 00548 output.attr("number",i+1); 00549 output.tag("NdMaterialOutput"); 00550 output.attr("classType", theMaterial[i]->getClassTag()); 00551 output.attr("tag", theMaterial[i]->getTag()); 00552 00553 output.tag("ResponseType","sigma11"); 00554 output.tag("ResponseType","sigma22"); 00555 output.tag("ResponseType","sigma33"); 00556 output.tag("ResponseType","sigma23"); 00557 output.tag("ResponseType","sigma31"); 00558 output.tag("ResponseType","sigma12"); 00559 00560 output.endTag(); // NdMaterialOutput 00561 output.endTag(); // GaussPoint 00562 } 00563 00564 theResponse = new ElementResponse(this, 5, Vector(Num_TotalGaussPts*6) ); 00565 00566 } else if (strcmp(argv[0],"gausspoint") == 0 || strcmp(argv[0],"GaussPoint") == 0) 00567 theResponse = new ElementResponse(this, 6, Vector(Num_TotalGaussPts*Num_Dim) ); 00568 00569 output.endTag(); // ElementOutput 00570 00571 return theResponse; 00572 00573 } 00574 00575 //============================================================================= 00576 int EightNodeBrick_u_p_U::getResponse(int responseID, Information &eleInfo) 00577 { 00578 if (responseID == 1) 00579 return eleInfo.setVector(this->getResistingForce()); 00580 00581 else if (responseID == 5) { 00582 static Vector stresses(Num_TotalGaussPts*6); 00583 stresstensor sigma; 00584 int cnt = 0; 00585 int i; 00586 for (i=0; i<Num_TotalGaussPts; i++) { 00587 sigma = theMaterial[i]->getStressTensor(); 00588 stresses(cnt++) = sigma.cval(1,1); //xx 00589 stresses(cnt++) = sigma.cval(2,2); //yy 00590 stresses(cnt++) = sigma.cval(3,3); //zz 00591 stresses(cnt++) = sigma.cval(2,3); //yz 00592 stresses(cnt++) = sigma.cval(3,1); //zx 00593 stresses(cnt++) = sigma.cval(1,2); //xy 00594 } 00595 return eleInfo.setVector(stresses); 00596 } 00597 00598 else if (responseID == 6) { 00599 static Vector Gpts(Num_TotalGaussPts*Num_Dim); 00600 tensor GCoord; 00601 int cnt = 0; 00602 int i,j; 00603 GCoord = getGaussPts(); 00604 for (i=0; i<Num_TotalGaussPts; i++) { 00605 for (j=0; j<Num_Dim; j++) { 00606 Gpts(cnt++) = GCoord.cval(i+1,j+1); 00607 } 00608 } 00609 return eleInfo.setVector(Gpts); 00610 } 00611 00612 else if (responseID == 7) { 00613 static Vector Gpts(Num_TotalGaussPts*2); 00614 stresstensor sigma; 00615 int i; 00616 for (i=0; i<Num_TotalGaussPts; i++) { 00617 sigma = theMaterial[i]->getStressTensor(); 00618 Gpts(i*2 ) = sigma.p_hydrostatic(); 00619 Gpts(i*2 +1) = sigma.q_deviatoric(); 00620 } 00621 return eleInfo.setVector(Gpts); 00622 } 00623 00624 else 00625 return (-1); 00626 } 00627 00628 //============================================================================= 00629 void EightNodeBrick_u_p_U::Print(OPS_Stream &s, int flag) 00630 { 00631 s << "EightNodeBrick_u_p_U, element id: " << this->getTag() << "\n"; 00632 s << "Connected external nodes: " << connectedExternalNodes << "\n"; 00633 00634 s << "Node 1: " << connectedExternalNodes(0) << "\n"; 00635 s << "Node 2: " << connectedExternalNodes(1) << "\n"; 00636 s << "Node 3: " << connectedExternalNodes(2) << "\n"; 00637 s << "Node 4: " << connectedExternalNodes(3) << "\n"; 00638 s << "Node 5: " << connectedExternalNodes(4) << "\n"; 00639 s << "Node 6: " << connectedExternalNodes(5) << "\n"; 00640 s << "Node 7: " << connectedExternalNodes(6) << "\n"; 00641 s << "Node 8: " << connectedExternalNodes(7) << "\n"; 00642 00643 s << "Material model: " << "\n"; 00644 00645 int GP_c_r, GP_c_s, GP_c_t, where; 00646 00647 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts ; GP_c_r++ ) { 00648 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts ; GP_c_s++ ) { 00649 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts ; GP_c_t++ ) { 00650 where = (GP_c_r*Num_IntegrationPts+GP_c_s)*Num_IntegrationPts+GP_c_t; 00651 s << "\n where = " << where+1 << "\n"; 00652 s << " r= " << GP_c_r << " s= " << GP_c_s << " t= " << GP_c_t << "\n"; 00653 theMaterial[where]->Print(s); 00654 } 00655 } 00656 } 00657 00658 } 00659 00660 //====================================================================== 00661 int EightNodeBrick_u_p_U::update() 00662 { 00663 int ret = 0; 00664 00665 double r = 0.0; 00666 double s = 0.0; 00667 double t = 0.0; 00668 00669 int Tdisp_dim[] = {Num_Nodes, Num_Dof}; 00670 tensor total_displacements(2, Tdisp_dim, 0.0); 00671 int tdisp_dim[] = {Num_Nodes, Num_Dim}; 00672 tensor total_disp(2, tdisp_dim, 0.0); 00673 00674 int dh_dim[] = {Num_Nodes, Num_Dim}; 00675 tensor dh(2, dh_dim, 0.0); 00676 00677 straintensor eps; 00678 00679 tensor dhGlobal; 00680 00681 total_displacements = getNodesDisp(); 00682 int i; 00683 for (i=1; i<=Num_Nodes; i++) { 00684 total_disp.val(i,1) = total_displacements.cval(i,1); 00685 total_disp.val(i,2) = total_displacements.cval(i,2); 00686 total_disp.val(i,3) = total_displacements.cval(i,3); 00687 } 00688 00689 int GP_c_r, GP_c_s, GP_c_t, where; 00690 00691 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 00692 r = pts[GP_c_r]; 00693 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 00694 s = pts[GP_c_s]; 00695 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 00696 t = pts[GP_c_t]; 00697 where = (GP_c_r*Num_IntegrationPts+GP_c_s)*Num_IntegrationPts+GP_c_t; 00698 dhGlobal = dh_Global(r,s,t); 00699 eps = total_disp("ia") * dhGlobal("ib"); 00700 eps.null_indices(); 00701 eps.symmetrize11(); 00702 if ( (theMaterial[where]->setTrialStrain(eps) ) ) 00703 opserr << "TwentyNodeBrick_u_p_U::update (tag: " << this->getTag() << "), not converged\n"; 00704 } 00705 } 00706 } 00707 00708 return ret; 00709 } 00710 00711 //====================================================================== 00712 const Vector& EightNodeBrick_u_p_U::getInternalForce () 00713 { 00714 static Vector Pforce(Num_ElemDof); 00715 00716 double r = 0.0; 00717 double rw = 0.0; 00718 double s = 0.0; 00719 double sw = 0.0; 00720 double t = 0.0; 00721 double tw = 0.0; 00722 double weight = 0.0; 00723 double det_of_Jacobian = 0.0; 00724 int where = 0; 00725 00726 int dh_dim[] = {Num_Nodes,Num_Dim}; 00727 tensor dh(2, dh_dim, 0.0); 00728 tensor Pins(2, dh_dim, 0.0); 00729 00730 tensor Jacobian; 00731 stresstensor Mstress; 00732 00733 int GP_c_r, GP_c_s, GP_c_t; 00734 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 00735 r = pts[GP_c_r]; 00736 rw = wts[GP_c_r]; 00737 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 00738 s = pts[GP_c_s]; 00739 sw = wts[GP_c_s]; 00740 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 00741 t = pts[GP_c_t]; 00742 tw = wts[GP_c_t]; 00743 where = (GP_c_r*Num_IntegrationPts+GP_c_s)*Num_IntegrationPts+GP_c_t; 00744 Jacobian = Jacobian_3D(r,s,t); 00745 det_of_Jacobian = Jacobian.determinant(); 00746 dh = dh_Global(r,s,t); 00747 weight = rw * sw * tw * det_of_Jacobian; 00748 Mstress = this->theMaterial[where]->getStressTensor(); 00749 Pins += (dh("Kj")*Mstress("ij"))*weight; 00750 } 00751 } 00752 } 00753 00754 Pforce.Zero(); // necessary 00755 00756 int i, j; 00757 for (i=0; i<Num_Nodes; i++) { 00758 for (j=0; j<Num_Dim; j++) { 00759 Pforce(i*Num_Dof+j) = Pins.cval(i+1, j+1); 00760 } 00761 } 00762 00763 return Pforce; 00764 } 00765 00766 00767 //====================================================================== 00768 const Vector& EightNodeBrick_u_p_U::getBodyForce () 00769 { 00770 static Vector Pforce(Num_ElemDof); 00771 00772 double r = 0.0; 00773 double rw = 0.0; 00774 double s = 0.0; 00775 double sw = 0.0; 00776 double t = 0.0; 00777 double tw = 0.0; 00778 double weight = 0.0; 00779 double det_of_Jacobian = 0.0; 00780 00781 int hp_dim[] = {Num_Nodes}; 00782 tensor hp(1, hp_dim, 0.0); 00783 tensor Pexf(1, hp_dim, 0.0); 00784 00785 tensor Jacobian; 00786 00787 int GP_c_r, GP_c_s, GP_c_t; 00788 00789 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 00790 r = pts[GP_c_r]; 00791 rw = wts[GP_c_r]; 00792 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 00793 s = pts[GP_c_s]; 00794 sw = wts[GP_c_s]; 00795 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 00796 t = pts[GP_c_t]; 00797 tw = wts[GP_c_t]; 00798 hp = shapeFunction(r,s,t); 00799 Jacobian = Jacobian_3D(r,s,t); 00800 det_of_Jacobian = Jacobian.determinant(); 00801 weight = rw * sw * tw * det_of_Jacobian; 00802 Pexf += hp *weight; 00803 } 00804 } 00805 } 00806 00807 Pforce.Zero(); // necessary 00808 00809 int i, j; 00810 for (i=0; i<Num_Nodes; i++) { 00811 for (j=0; j<Num_Dim; j++) { 00812 Pforce(i*Num_Dof +j) = Pexf.cval(i+1) * bf(j) * (1.0-poro) * rho_s; 00813 Pforce(i*Num_Dof +j +4) = Pexf.cval(i+1) * bf(j) * poro * rho_f; 00814 } 00815 } 00816 00817 return Pforce; 00818 } 00819 00820 00821 //====================================================================== 00822 const Matrix& EightNodeBrick_u_p_U::getStiff00 (void) 00823 { 00824 MCK.Zero(); // necessary 00825 00826 int i, j, m; 00827 00828 tensor tG = getStiffnessTensorG12(); 00829 00830 //G1 and G1^T, Note *(alpha-poro) here! 00831 for ( i=0 ; i<Num_Nodes; i++ ) { 00832 for ( j=0; j<Num_Nodes; j++ ) { 00833 for( m=0; m<Num_Dim; m++) { 00834 MCK(i*Num_Dof+m, j*Num_Dof+3) = -tG.cval(i+1, m+1, j+1) *(alpha-poro); 00835 MCK(j*Num_Dof+3, i*Num_Dof+m) = -tG.cval(i+1, m+1, j+1) *(alpha-poro); 00836 } 00837 } 00838 } 00839 00840 //G2 and G2^T, Note *poro here! 00841 for ( i=0 ; i<Num_Nodes; i++ ) { 00842 for ( j=0; j<Num_Nodes; j++ ) { 00843 for( m=0; m<Num_Dim; m++) { 00844 MCK(i*Num_Dof+m+4, j*Num_Dof+3) = -tG.cval(i+1, m+1, j+1) *poro; 00845 MCK(j*Num_Dof+3, i*Num_Dof+m+4) = -tG.cval(i+1, m+1, j+1) *poro; 00846 } 00847 } 00848 } 00849 00850 00851 //P 00852 if (ks == 0.0 || kf == 0.0) { 00853 opserr<<" Error, EightNodeBrick_u_p_U::getStiffnessTensorP -- solid and/or fluid bulk modulus is zero\n"; 00854 exit(-1); 00855 } 00856 00857 double oneOverQ = poro/kf + (alpha-poro)/ks; 00858 00859 tensor tP = getMassTensorMsf(); 00860 00861 for ( i=0 ; i<Num_Nodes; i++ ) { 00862 for ( j=0; j<Num_Nodes; j++ ) { 00863 MCK(i*Num_Dof+3, j*Num_Dof+3) = tP.cval(i+1, j+1) * (-oneOverQ); 00864 } 00865 } 00866 00867 return MCK; 00868 } 00869 00870 //====================================================================== 00871 const Matrix& EightNodeBrick_u_p_U::getStiff (int Ki_flag) 00872 { 00873 if (Ki_flag != 0 && Ki_flag != 1) { 00874 opserr << "Error EightNodeBrick_u_p_U::getStiff() - illegal use\n"; 00875 exit(-1); 00876 } 00877 00878 if (Ki_flag == 0 && Ki != 0) 00879 return *Ki; 00880 00881 tensor tKep = getStiffnessTensorKep(); 00882 00883 int i, j, m, n; 00884 00885 this->getStiff00(); 00886 00887 // + Kep 00888 for ( i=0 ; i<Num_Nodes; i++ ) { 00889 for ( j=0; j<Num_Nodes; j++ ) { 00890 for( m=0; m<Num_Dim; m++) { 00891 for( n=0; n<Num_Dim; n++) 00892 MCK(i*Num_Dof+m, j*Num_Dof+n) = tKep.cval(i+1, m+1, n+1, j+1); 00893 } 00894 } 00895 } 00896 00897 if( Ki_flag == 1) 00898 return MCK; 00899 00900 Ki = new Matrix(MCK); 00901 00902 if (Ki == 0) { 00903 opserr << "Error EightNodeBrick_u_p_U::getStiff() -"; 00904 opserr << "ran out of memory\n"; 00905 exit(-1); 00906 } 00907 00908 return *Ki; 00909 } 00910 00911 //====================================================================== 00912 tensor EightNodeBrick_u_p_U::getStiffnessTensorKep( ) 00913 { 00914 int K_dim[] = {Num_Nodes, Num_Dim, Num_Dim, Num_Nodes}; 00915 tensor Kep(4, K_dim, 0.0); 00916 tensor Kkt(4, K_dim, 0.0); 00917 00918 double r = 0.0; 00919 double rw = 0.0; 00920 double s = 0.0; 00921 double sw = 0.0; 00922 double t = 0.0; 00923 double tw = 0.0; 00924 int where = 0; 00925 double weight = 0.0; 00926 double det_of_Jacobian = 0.0; 00927 00928 int dh_dim[] = {Num_Nodes,Num_Dim}; 00929 tensor dh(2, dh_dim, 0.0); 00930 00931 tensor Constitutive; 00932 00933 tensor Jacobian; 00934 tensor dhGlobal; 00935 00936 int GP_c_r, GP_c_s, GP_c_t; 00937 00938 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 00939 r = pts[GP_c_r]; 00940 rw = wts[GP_c_r]; 00941 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 00942 s = pts[GP_c_s]; 00943 sw = wts[GP_c_s]; 00944 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 00945 t = pts[GP_c_t]; 00946 tw = wts[GP_c_t]; 00947 where = (GP_c_r*Num_IntegrationPts+GP_c_s)*Num_IntegrationPts+GP_c_t; 00948 Jacobian = Jacobian_3D(r,s,t); 00949 det_of_Jacobian = Jacobian.determinant(); 00950 dhGlobal = dh_Global(r,s,t); 00951 weight = rw * sw * tw * det_of_Jacobian; 00952 Constitutive = theMaterial[where]->getTangentTensor(); 00953 Kkt = dhGlobal("kj")*Constitutive("ijml"); 00954 Kkt = Kkt("kiml")*dhGlobal("pl")*weight; 00955 Kep = Kep + Kkt; 00956 } 00957 } 00958 } 00959 00960 return Kep; 00961 } 00962 00963 //====================================================================== 00964 tensor EightNodeBrick_u_p_U::getStiffnessTensorG12() 00965 { 00966 // This is for G1 and G2 00967 // G1 = (alpha-poro) *G; 00968 // G2 = poro *G; 00969 00970 int G_dim[] = {Num_Nodes, Num_Dim, Num_Nodes}; 00971 tensor G(3, G_dim, 0.0); 00972 00973 double r = 0.0; 00974 double rw = 0.0; 00975 double s = 0.0; 00976 double sw = 0.0; 00977 double t = 0.0; 00978 double tw = 0.0; 00979 double weight = 0.0; 00980 00981 int dh_dim[] = {Num_Nodes,Num_Dim}; 00982 tensor dh(2, dh_dim, 0.0); 00983 int hp_dim[] = {Num_Nodes}; 00984 tensor hp(1, hp_dim, 0.0); 00985 00986 double det_of_Jacobian = 0.0; 00987 00988 tensor Jacobian; 00989 tensor dhGlobal; 00990 00991 int GP_c_r, GP_c_s, GP_c_t; 00992 00993 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 00994 r = pts[GP_c_r]; 00995 rw = wts[GP_c_r]; 00996 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 00997 s = pts[GP_c_s]; 00998 sw = wts[GP_c_s]; 00999 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 01000 t = pts[GP_c_t]; 01001 tw = wts[GP_c_t]; 01002 hp= shapeFunction(r,s,t); 01003 Jacobian = Jacobian_3D(r,s,t); 01004 dhGlobal = dh_Global(r,s,t); 01005 det_of_Jacobian = Jacobian.determinant(); 01006 weight = rw * sw * tw * det_of_Jacobian; 01007 G += dhGlobal("ki")*hp("m") * weight; 01008 } 01009 } 01010 } 01011 01012 return G; 01013 } 01014 01015 //====================================================================== 01016 tensor EightNodeBrick_u_p_U::getDampTensorC123() 01017 { 01018 // This is for C1, C2 and C3, C1 = C2 = c3 01019 // Since solid and fluid shape function the same 01020 01021 int perm_dim[] = {Num_Dim, Num_Dim}; 01022 tensor perm_inv(2, perm_dim, 0.0); 01023 01024 perm_inv.val(1,1) = 1.0/perm(0); 01025 perm_inv.val(2,2) = 1.0/perm(1); 01026 perm_inv.val(3,3) = 1.0/perm(2); 01027 01028 int C_dim[] = {Num_Nodes,Num_Dim,Num_Dim,Num_Nodes}; 01029 tensor C123(4,C_dim,0.0); 01030 int c_dim[] = {Num_Nodes,Num_Dim,Num_Dim}; 01031 tensor c123(3,c_dim,0.0); 01032 01033 double r = 0.0; 01034 double rw = 0.0; 01035 double s = 0.0; 01036 double sw = 0.0; 01037 double t = 0.0; 01038 double tw = 0.0; 01039 double weight = 0.0; 01040 double det_of_Jacobian = 0.0; 01041 01042 int hp_dim[] = {Num_Nodes}; 01043 tensor hp(1, hp_dim,0.0); 01044 01045 tensor Jacobian; 01046 01047 int GP_c_r, GP_c_s, GP_c_t; 01048 01049 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 01050 r = pts[GP_c_r]; 01051 rw = wts[GP_c_r]; 01052 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 01053 s = pts[GP_c_s]; 01054 sw = wts[GP_c_s]; 01055 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 01056 t = pts[GP_c_t]; 01057 tw = wts[GP_c_t]; 01058 hp = shapeFunction(r,s,t); 01059 Jacobian = Jacobian_3D(r,s,t); 01060 det_of_Jacobian = Jacobian.determinant(); 01061 weight = rw * sw * tw * det_of_Jacobian; 01062 c123 = hp("k")*perm_inv("ij"); 01063 C123 += c123("kij")*hp("m") *weight; 01064 } 01065 } 01066 } 01067 01068 return C123; 01069 } 01070 01071 //====================================================================== 01072 tensor EightNodeBrick_u_p_U::getMassTensorMsf() 01073 { 01074 // This is for Ms and Mf -> M_kl 01075 // Ms = Msf * (1.0-poro)*rho_s 01076 // Mf = Msf * poro*rho_f 01077 01078 // Also this is for Compression term Pc 01079 // Pc = Msf * oneOverQ 01080 // = Msf * (poro/kf + (alpha-poro)/ks) 01081 01082 int M_dim[] = {Num_Nodes, Num_Nodes}; 01083 tensor Msf(2, M_dim, 0.0); 01084 01085 double r = 0.0; 01086 double rw = 0.0; 01087 double s = 0.0; 01088 double sw = 0.0; 01089 double t = 0.0; 01090 double tw = 0.0; 01091 double weight = 0.0; 01092 double det_of_Jacobian = 0.0; 01093 01094 int dh_dim[] = {Num_Nodes,Num_Dim}; 01095 tensor dh(2, dh_dim, 0.0); 01096 int hp_dim[] = {Num_Nodes}; 01097 tensor hp(1, hp_dim,0.0); 01098 01099 tensor Jacobian; 01100 01101 int GP_c_r, GP_c_s, GP_c_t; 01102 01103 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 01104 r = pts[GP_c_r]; 01105 rw = wts[GP_c_r]; 01106 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 01107 s = pts[GP_c_s]; 01108 sw = wts[GP_c_s]; 01109 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 01110 t = pts[GP_c_t]; 01111 tw = wts[GP_c_t]; 01112 hp = shapeFunction(r,s,t); 01113 Jacobian = Jacobian_3D(r,s,t); 01114 det_of_Jacobian = Jacobian.determinant(); 01115 weight = rw * sw * tw * det_of_Jacobian; 01116 Msf += hp("m")*hp("n")*weight; 01117 } 01118 } 01119 } 01120 01121 return Msf; 01122 } 01123 01124 //====================================================================== 01125 tensor EightNodeBrick_u_p_U::Jacobian_3D(double x, double y, double z) 01126 { 01127 tensor N_C = this->getNodesCrds(); 01128 tensor dh = this->shapeFunctionDerivative(x, y, z); 01129 01130 tensor J3D = N_C("ki") * dh("kj"); 01131 J3D.null_indices(); 01132 return J3D; 01133 } 01134 01135 //====================================================================== 01136 tensor EightNodeBrick_u_p_U::Jacobian_3Dinv(double x, double y, double z) 01137 { 01138 tensor J = this->Jacobian_3D(x,y,z); 01139 return J.inverse(); 01140 } 01141 01142 //====================================================================== 01143 tensor EightNodeBrick_u_p_U::dh_Global(double x, double y, double z) 01144 { 01145 tensor JacobianINV0 = this->Jacobian_3Dinv(x, y, z); 01146 tensor dh = this->shapeFunctionDerivative(x, y, z); 01147 tensor dhGlobal = dh("ik") * JacobianINV0("kj"); 01148 dhGlobal.null_indices(); 01149 return dhGlobal; 01150 } 01151 01152 //====================================================================== 01153 tensor EightNodeBrick_u_p_U::getNodesCrds(void) 01154 { 01155 int i,j; 01156 int dimX[] = {Num_Nodes,Num_Dim}; 01157 tensor N_coord(2, dimX, 0.0); 01158 01159 for (i=0; i<Num_Nodes; i++) { 01160 const Vector&TNodesCrds = theNodes[i]->getCrds(); 01161 for (j=0; j<Num_Dim; j++) { 01162 N_coord.val(i+1,j+1) = TNodesCrds(j); 01163 } 01164 } 01165 01166 return N_coord; 01167 } 01168 01169 //====================================================================== 01170 tensor EightNodeBrick_u_p_U::getNodesDisp(void) 01171 { 01172 int i,j; 01173 int dimU[] = {Num_Nodes, Num_Dof}; 01174 tensor total_disp(2, dimU, 0.0); 01175 01176 for (i=0; i<Num_Nodes; i++) { 01177 const Vector&TNodesDisp = theNodes[i]->getTrialDisp(); 01178 for (j=0; j<Num_Dof; j++) { 01179 total_disp.val(i+1,j+1) = TNodesDisp(j); 01180 } 01181 } 01182 01183 return total_disp; 01184 } 01185 01186 //====================================================================== 01187 double EightNodeBrick_u_p_U::getPorePressure(double x1, double x2, double x3) 01188 { 01189 double pp = 0.0; 01190 int i; 01191 01192 for (i=0; i<Num_Nodes; i++) { 01193 const Vector& T_disp = theNodes[i]->getTrialDisp(); 01194 pp += shapeFunction(x1,x2,x3).cval(i+1) * T_disp(3); 01195 } 01196 01197 return pp; 01198 } 01199 01200 //====================================================================== 01201 tensor EightNodeBrick_u_p_U::shapeFunction(double r1, double r2, double r3) 01202 { 01203 int Hfun[] = {Num_Nodes}; 01204 tensor h(1, Hfun, 0.0); 01205 01206 h.val(8)=(1.0+r1)*(1.0-r2)*(1.0-r3)*0.125; 01207 h.val(7)=(1.0-r1)*(1.0-r2)*(1.0-r3)*0.125; 01208 h.val(6)=(1.0-r1)*(1.0+r2)*(1.0-r3)*0.125; 01209 h.val(5)=(1.0+r1)*(1.0+r2)*(1.0-r3)*0.125; 01210 h.val(4)=(1.0+r1)*(1.0-r2)*(1.0+r3)*0.125; 01211 h.val(3)=(1.0-r1)*(1.0-r2)*(1.0+r3)*0.125; 01212 h.val(2)=(1.0-r1)*(1.0+r2)*(1.0+r3)*0.125; 01213 h.val(1)=(1.0+r1)*(1.0+r2)*(1.0+r3)*0.125; 01214 01215 return h; 01216 } 01217 01218 01219 //============================================================== 01220 tensor EightNodeBrick_u_p_U::shapeFunctionDerivative(double r1, double r2, double r3) 01221 { 01222 int DHfun[] = {Num_Nodes, Num_Dim}; 01223 tensor dh(2, DHfun, 0.0); 01224 01225 // node number 8 01226 dh.val(8,1)= (1.0-r2)*(1.0-r3)*0.125; 01227 dh.val(8,2)=-(1.0+r1)*(1.0-r3)*0.125; 01228 dh.val(8,3)=-(1.0+r1)*(1.0-r2)*0.125; 01229 // node number 7 01230 dh.val(7,1)=-(1.0-r2)*(1.0-r3)*0.125; 01231 dh.val(7,2)=-(1.0-r1)*(1.0-r3)*0.125; 01232 dh.val(7,3)=-(1.0-r1)*(1.0-r2)*0.125; 01233 // node number 6 01234 dh.val(6,1)=-(1.0+r2)*(1.0-r3)*0.125; 01235 dh.val(6,2)= (1.0-r1)*(1.0-r3)*0.125; 01236 dh.val(6,3)=-(1.0-r1)*(1.0+r2)*0.125; 01237 // node number 5 01238 dh.val(5,1)= (1.0+r2)*(1.0-r3)*0.125; 01239 dh.val(5,2)= (1.0+r1)*(1.0-r3)*0.125; 01240 dh.val(5,3)=-(1.0+r1)*(1.0+r2)*0.125; 01241 // node number 4 01242 dh.val(4,1)= (1.0-r2)*(1.0+r3)*0.125; 01243 dh.val(4,2)=-(1.0+r1)*(1.0+r3)*0.125; 01244 dh.val(4,3)= (1.0+r1)*(1.0-r2)*0.125; 01245 // node number 3 01246 dh.val(3,1)=-(1.0-r2)*(1.0+r3)*0.125; 01247 dh.val(3,2)=-(1.0-r1)*(1.0+r3)*0.125; 01248 dh.val(3,3)= (1.0-r1)*(1.0-r2)*0.125; 01249 // node number 2 01250 dh.val(2,1)=-(1.0+r2)*(1.0+r3)*0.125; 01251 dh.val(2,2)= (1.0-r1)*(1.0+r3)*0.125; 01252 dh.val(2,3)= (1.0-r1)*(1.0+r2)*0.125; 01253 // node number 1 01254 dh.val(1,1)= (1.0+r2)*(1.0+r3)*0.125; 01255 dh.val(1,2)= (1.0+r1)*(1.0+r3)*0.125; 01256 dh.val(1,3)= (1.0+r1)*(1.0+r2)*0.125; 01257 01258 return dh; 01259 } 01260 01261 //============================================================== 01262 tensor EightNodeBrick_u_p_U::getGaussPts(void) 01263 { 01264 int dimensions1[] = {Num_TotalGaussPts, Num_Dim}; 01265 tensor Gs(2, dimensions1, 0.0); 01266 int dimensions2[] = {Num_Nodes}; 01267 tensor shp(1, dimensions2, 0.0); 01268 01269 double r = 0.0; 01270 double s = 0.0; 01271 double t = 0.0; 01272 int i, j, where; 01273 01274 int GP_c_r, GP_c_s, GP_c_t; 01275 01276 for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) { 01277 r = pts[GP_c_r]; 01278 for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) { 01279 s = pts[GP_c_s]; 01280 for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) { 01281 t = pts[GP_c_t]; 01282 where = (GP_c_r*Num_IntegrationPts+GP_c_s)*Num_IntegrationPts+GP_c_t; 01283 shp = shapeFunction(r,s,t); 01284 for (i=0; i<Num_Nodes; i++) { 01285 const Vector& T_Crds = theNodes[i]->getCrds(); 01286 for (j=0; j<Num_Dim; j++) { 01287 Gs.val(where+1, j+1) += shp.cval(i+1) * T_Crds(j); 01288 } 01289 } 01290 } 01291 } 01292 } 01293 01294 return Gs; 01295 } 01296 01297 01298 #endif 01299 01300 01301
Generated on Mon Oct 23 15:05:12 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents