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