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