OpenSees » SRC » element » UP_ucdavis

EightNode_LDBrick_u_p.cpp

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

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