My title

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.3 $
00022 // $Date: 2001/05/22 07:33:23 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/material/section/FiberSection3d.cpp,v $
00024                                                                         
00025 // Written: fmk
00026 // Created: 04/04
00027 //
00028 // Description: This file contains the class implementation of FiberSection2d.
00029 
00030 #include <stdlib.h>
00031 
00032 #include <Channel.h>
00033 #include <Vector.h>
00034 #include <Matrix.h>
00035 #include <MatrixUtil.h>
00036 #include <Fiber.h>
00037 #include <classTags.h>
00038 #include <FiberSection3d.h>
00039 #include <ID.h>
00040 #include <FEM_ObjectBroker.h>
00041 #include <Information.h>
00042 #include <MaterialResponse.h>
00043 #include <UniaxialMaterial.h>
00044 
00045 ID FiberSection3d::code(3);
00046 
00047 // constructors:
00048 FiberSection3d::FiberSection3d(int tag, int num, Fiber **fibers): 
00049   SectionForceDeformation(tag, SEC_TAG_FiberSection3d),
00050   numFibers(num), theMaterials(0), matData(0), e(3), eCommit(3), s(0), ks(0)
00051 {
00052   if (numFibers != 0) {
00053     theMaterials = new UniaxialMaterial *[numFibers];
00054 
00055     if (theMaterials == 0)
00056       g3ErrorHandler->fatal("%s -- failed to allocate Material pointers",
00057        "FiberSection3d::FiberSection3d");
00058 
00059     matData = new double [numFibers*3];
00060 
00061     if (matData == 0)
00062       g3ErrorHandler->fatal("%s -- failed to allocate double array for material data",
00063        "FiberSection3d::FiberSection3d");
00064     
00065     for (int i = 0; i < numFibers; i++) {
00066       Fiber *theFiber = fibers[i];
00067       double yLoc, zLoc, Area;
00068       theFiber->getFiberLocation(yLoc, zLoc);
00069       Area = theFiber->getArea();
00070 
00071       matData[i*3] = -yLoc;
00072       matData[i*3+1] = zLoc;
00073       matData[i*3+2] = Area;
00074       UniaxialMaterial *theMat = theFiber->getMaterial();
00075       theMaterials[i] = theMat->getCopy();
00076 
00077       if (theMaterials[i] == 0)
00078  g3ErrorHandler->fatal("%s -- failed to get copy of a Material",
00079          "FiberSection3d::FiberSection3d");
00080     }    
00081   }
00082 
00083   s = new Vector(sData, 3);
00084   ks = new Matrix(kData, 3, 3);
00085 
00086   sData[0] = 0.0;
00087   sData[1] = 0.0;
00088   sData[2] = 0.0;
00089 
00090   for (int i=0; i<9; i++)
00091     kData[i] = 0.0;
00092 
00093   code(0) = SECTION_RESPONSE_P;
00094   code(1) = SECTION_RESPONSE_MZ;
00095   code(2) = SECTION_RESPONSE_MY;
00096 }
00097 
00098 // constructor for blank object that recvSelf needs to be invoked upon
00099 FiberSection3d::FiberSection3d():
00100   SectionForceDeformation(0, SEC_TAG_FiberSection3d),
00101   numFibers(0), theMaterials(0), matData(0), e(3), eCommit(3), s(0), ks(0)
00102 {
00103   s = new Vector(sData, 3);
00104   ks = new Matrix(kData, 3, 3);
00105 
00106   sData[0] = 0.0;
00107   sData[1] = 0.0;
00108   sData[2] = 0.0;
00109 
00110   for (int i=0; i<9; i++)
00111     kData[i] = 0.0;
00112 
00113   code(0) = SECTION_RESPONSE_P;
00114   code(1) = SECTION_RESPONSE_MZ;
00115   code(2) = SECTION_RESPONSE_MY;
00116 
00117 }
00118 
00119 int
00120 FiberSection3d::addFiber(Fiber &newFiber)
00121 {
00122   // need to create a larger array
00123   int newSize = numFibers+1;
00124 
00125   UniaxialMaterial **newArray = new UniaxialMaterial *[newSize]; 
00126   double *newMatData = new double [3 * newSize];
00127   
00128   if (newArray == 0 || newMatData == 0) {
00129     g3ErrorHandler->fatal("%s -- failed to allocate Fiber pointers",
00130      "FiberSection3d::addFiber");
00131     return -1;
00132   }
00133 
00134   // copy the old pointers
00135   for (int i = 0; i < numFibers; i++) {
00136     newArray[i] = theMaterials[i];
00137     newMatData[3*i] = matData[3*i];
00138     newMatData[3*i+1] = matData[3*i+1];
00139     newMatData[3*i+2] = matData[3*i+2];
00140   }
00141   // set the new pointers
00142   double yLoc, zLoc, Area;
00143   newFiber.getFiberLocation(yLoc, zLoc);
00144   Area = newFiber.getArea();
00145   newMatData[numFibers*3] = -yLoc;
00146   newMatData[numFibers*3+1] = zLoc;
00147   newMatData[numFibers*3+2] = Area;
00148   UniaxialMaterial *theMat = newFiber.getMaterial();
00149   newArray[numFibers] = theMat->getCopy();
00150 
00151   if (newArray[numFibers] == 0) {
00152     g3ErrorHandler->fatal("%s -- failed to get copy of a Material",
00153      "FiberSection3d::addFiber");
00154 
00155     delete [] newArray;
00156     delete [] newMatData;
00157     return -1;
00158   }
00159 
00160   numFibers++;
00161   
00162   if (theMaterials != 0) {
00163     delete [] theMaterials;
00164     delete [] matData;
00165   }
00166 
00167   theMaterials = newArray;
00168   matData = newMatData;
00169 
00170   return 0;
00171 }
00172 
00173 
00174 
00175 // destructor:
00176 FiberSection3d::~FiberSection3d()
00177 {
00178   if (theMaterials != 0) {
00179     for (int i = 0; i < numFibers; i++)
00180       if (theMaterials[i] != 0)
00181  delete theMaterials[i];
00182       
00183     delete [] theMaterials;
00184   }
00185 
00186   if (matData != 0)
00187     delete [] matData;
00188 
00189   if (s != 0)
00190     delete s;
00191 
00192   if (ks != 0)
00193     delete ks;
00194 }
00195 
00196 int FiberSection3d::setTrialSectionDeformation (const Vector &deforms)
00197 {
00198   int res = 0;
00199   e = deforms;
00200 
00201   kData[0] = 0.0; kData[1] = 0.0; kData[2] = 0.0; kData[3] = 0.0;
00202   kData[4] = 0.0; kData[5] = 0.0; kData[6] = 0.0; kData[7] = 0.0;
00203   kData[8] = 0.0; 
00204   sData[0] = 0.0; sData[1] = 0.0;  sData[2] = 0.0; 
00205 
00206   int loc = 0;
00207 
00208   double d0 = deforms(0);
00209   double d1 = deforms(1);
00210   double d2 = deforms(2);
00211 
00212   for (int i = 0; i < numFibers; i++) {
00213     UniaxialMaterial *theMat = theMaterials[i];
00214     double y = matData[loc++];
00215     double z = matData[loc++];
00216     double A = matData[loc++];
00217 
00218     // determine material strain and set it
00219     double strain = d0 + y*d1 + z*d2;
00220     double tangent, stress;
00221     res = theMat->setTrial(strain, stress, tangent);
00222 
00223     double value = tangent * A;
00224     double vas1 = y*value;
00225     double vas2 = z*value;
00226     double vas1as2 = vas1*z;
00227 
00228     kData[0] += value;
00229     kData[1] += vas1;
00230     kData[2] += vas2;
00231     
00232     kData[4] += vas1 * y;
00233     kData[5] += vas1as2;
00234     
00235     kData[8] += vas2 * z; 
00236 
00237     double fs0 = stress * A;
00238     sData[0] += fs0;
00239     sData[1] += fs0 * y;
00240     sData[2] += fs0 * z;
00241   }
00242 
00243   kData[3] = kData[1];
00244   kData[6] = kData[2];
00245   kData[7] = kData[5];
00246 
00247   return res;
00248 }
00249 
00250 const Vector&
00251 FiberSection3d::getSectionDeformation(void)
00252 {
00253   return e;
00254 }
00255 
00256 const Matrix&
00257 FiberSection3d::getSectionTangent(void)
00258 {
00259   return *ks;
00260 }
00261 
00262 const Matrix&
00263 FiberSection3d::getSectionSecant(void)
00264 {
00265   return *ks;
00266 }
00267 
00268 const Vector&
00269 FiberSection3d::getStressResultant(void)
00270 {
00271   return *s;
00272 }
00273 
00274 SectionForceDeformation*
00275 FiberSection3d::getCopy(void)
00276 {
00277   FiberSection3d *theCopy = new FiberSection3d ();
00278   theCopy->setTag(this->getTag());
00279 
00280   theCopy->numFibers = numFibers;
00281 
00282   if (numFibers != 0) {
00283     theCopy->theMaterials = new UniaxialMaterial *[numFibers];
00284 
00285     if (theCopy->theMaterials == 0)
00286       g3ErrorHandler->fatal("%s -- failed to allocate Material pointers",
00287        "FiberSection3d::FiberSection3d");
00288 
00289     theCopy->matData = new double [numFibers*3];
00290 
00291     if (theCopy->matData == 0)
00292       g3ErrorHandler->fatal("%s -- failed to allocate double array for material data",
00293        "FiberSection3d::FiberSection3d");
00294     
00295     for (int i = 0; i < numFibers; i++) {
00296       theCopy->matData[i*3] = matData[i*3];
00297       theCopy->matData[i*3+1] = matData[i*3+1];
00298       theCopy->matData[i*3+2] = matData[i*3+2];
00299       theCopy->theMaterials[i] = theMaterials[i]->getCopy();
00300 
00301       if (theCopy->theMaterials[i] == 0)
00302  g3ErrorHandler->fatal("%s -- failed to get copy of a Material",
00303          "FiberSection3d::getCopy");
00304     }    
00305   }
00306 
00307   theCopy->eCommit = eCommit;
00308   theCopy->e = e;
00309 
00310   for (int i=0; i<9; i++)
00311     theCopy->kData[i] = kData[i];
00312 
00313   theCopy->sData[0] = sData[0];
00314   theCopy->sData[1] = sData[1];
00315   theCopy->sData[2] = sData[2];
00316 
00317   theCopy->s = new Vector(theCopy->sData, 3);
00318   theCopy->ks = new Matrix(theCopy->kData, 3, 3);
00319 
00320   return theCopy;
00321 }
00322 
00323 const ID&
00324 FiberSection3d::getType () const
00325 {
00326   return code;
00327 }
00328 
00329 int
00330 FiberSection3d::getOrder () const
00331 {
00332   return 3;
00333 }
00334 
00335 int
00336 FiberSection3d::commitState(void)
00337 {
00338   int err = 0;
00339 
00340   for (int i = 0; i < numFibers; i++)
00341     err += theMaterials[i]->commitState();
00342 
00343   eCommit = e;
00344 
00345   return err;
00346 }
00347 
00348 int
00349 FiberSection3d::revertToLastCommit(void)
00350 {
00351   int err = 0;
00352 
00353   // Last committed section deformations
00354   e = eCommit;
00355 
00356 
00357   kData[0] = 0.0; kData[1] = 0.0; kData[2] = 0.0; kData[3] = 0.0;
00358   kData[4] = 0.0; kData[5] = 0.0; kData[6] = 0.0; kData[7] = 0.0;
00359   kData[8] = 0.0; 
00360   sData[0] = 0.0; sData[1] = 0.0;  sData[2] = 0.0; 
00361 
00362   int loc = 0;
00363 
00364   for (int i = 0; i < numFibers; i++) {
00365     UniaxialMaterial *theMat = theMaterials[i];
00366     double y = matData[loc++];
00367     double z = matData[loc++];
00368     double A = matData[loc++];
00369 
00370     // invoke revertToLast on the material
00371     err += theMat->revertToLastCommit();
00372 
00373     double tangent = theMat->getTangent();
00374     double stress = theMat->getStress();
00375 
00376     double value = tangent * A;
00377     double vas1 = y*value;
00378     double vas2 = z*value;
00379     double vas1as2 = vas1*z;
00380 
00381     kData[0] += value;
00382     kData[1] += vas1;
00383     kData[2] += vas2;
00384     
00385     kData[4] += vas1 * y;
00386     kData[5] += vas1as2;
00387     
00388     kData[8] += vas2 * z; 
00389 
00390     double fs0 = stress * A;
00391     sData[0] += fs0;
00392     sData[1] += fs0 * y;
00393     sData[2] += fs0 * z;
00394   }
00395 
00396   kData[3] = kData[1];
00397   kData[6] = kData[2];
00398   kData[7] = kData[5];
00399 
00400   return err;
00401 }
00402 
00403 int
00404 FiberSection3d::revertToStart(void)
00405 {
00406   // revert the fibers to start    
00407   int err = 0;
00408 
00409 
00410   kData[0] = 0.0; kData[1] = 0.0; kData[2] = 0.0; kData[3] = 0.0;
00411   kData[4] = 0.0; kData[5] = 0.0; kData[6] = 0.0; kData[7] = 0.0;
00412   kData[8] = 0.0; 
00413   sData[0] = 0.0; sData[1] = 0.0;  sData[2] = 0.0; 
00414 
00415   int loc = 0;
00416 
00417   for (int i = 0; i < numFibers; i++) {
00418     UniaxialMaterial *theMat = theMaterials[i];
00419     double y = matData[loc++];
00420     double z = matData[loc++];
00421     double A = matData[loc++];
00422 
00423     // invoke revertToStart on the material
00424     err += theMat->revertToStart();
00425 
00426     double tangent = theMat->getTangent();
00427     double stress = theMat->getStress();
00428 
00429     double value = tangent * A;
00430     double vas1 = y*value;
00431     double vas2 = z*value;
00432     double vas1as2 = vas1*z;
00433 
00434     kData[0] += value;
00435     kData[1] += vas1;
00436     kData[2] += vas2;
00437     
00438     kData[4] += vas1 * y;
00439     kData[5] += vas1as2;
00440     
00441     kData[8] += vas2 * z; 
00442 
00443     double fs0 = stress * A;
00444     sData[0] += fs0;
00445     sData[1] += fs0 * y;
00446     sData[2] += fs0 * z;
00447   }
00448 
00449   kData[3] = kData[1];
00450   kData[6] = kData[2];
00451   kData[7] = kData[5];
00452 
00453   return err;
00454 }
00455 
00456 int
00457 FiberSection3d::sendSelf(int commitTag, Channel &theChannel)
00458 {
00459   int res = 0;
00460 
00461   // create an id to send objects tag and numFibers, 
00462   //     size 3 so no conflict with matData below if just 1 fiber
00463   static ID data(3);
00464   data(0) = this->getTag();
00465   data(1) = numFibers;
00466   int dbTag = this->getDbTag();
00467   res += theChannel.sendID(dbTag, commitTag, data);
00468   if (res < 0) {
00469     g3ErrorHandler->warning("%s - failed to send ID data",
00470        "FiberSection2d::sendSelf");
00471     return res;
00472   }    
00473 
00474   if (numFibers != 0) {
00475     
00476     // create an id containingg classTag and dbTag for each material & send it
00477     ID materialData(2*numFibers);
00478     for (int i=0; i<numFibers; i++) {
00479       UniaxialMaterial *theMat = theMaterials[i];
00480       materialData(2*i) = theMat->getClassTag();
00481       int matDbTag = theMat->getDbTag();
00482       if (matDbTag == 0) {
00483  matDbTag = theChannel.getDbTag();
00484  if (matDbTag != 0)
00485    theMat->setDbTag(matDbTag);
00486       }
00487       materialData(2*i+1) = matDbTag;
00488     }    
00489     
00490     res += theChannel.sendID(dbTag, commitTag, materialData);
00491     if (res < 0) {
00492       g3ErrorHandler->warning("%s - failed to send material data",
00493          "FiberSection2d::sendSelf");
00494       return res;
00495     }    
00496 
00497     // send the fiber data, i.e. area and loc
00498     Vector fiberData(matData, 3*numFibers);
00499     res += theChannel.sendVector(dbTag, commitTag, fiberData);
00500     if (res < 0) {
00501       g3ErrorHandler->warning("%s - failed to send material data",
00502          "FiberSection2d::sendSelf");
00503       return res;
00504     }    
00505 
00506     // now invoke send(0 on all the materials
00507     for (int j=0; j<numFibers; j++)
00508       theMaterials[j]->sendSelf(commitTag, theChannel);
00509   }
00510 
00511   return res;
00512 }
00513 
00514 int
00515 FiberSection3d::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00516 {
00517   int res = 0;
00518 
00519   static ID data(3);
00520   
00521   int dbTag = this->getDbTag();
00522   res += theChannel.recvID(dbTag, commitTag, data);
00523   if (res < 0) {
00524     g3ErrorHandler->warning("%s - failed to recv ID data",
00525        "FiberSection2d::recvSelf");
00526     return res;
00527   }    
00528   this->setTag(data(0));
00529 
00530   // recv data about materials objects, classTag and dbTag
00531   if (data(1) != 0) {
00532     ID materialData(2*data(1));
00533     res += theChannel.recvID(dbTag, commitTag, materialData);
00534     if (res < 0) {
00535       g3ErrorHandler->warning("%s - failed to send material data",
00536          "FiberSection2d::recvSelf");
00537       return res;
00538     }    
00539 
00540     // if current arrays not of correct size, release old and resize
00541     if (theMaterials == 0 || numFibers != data(1)) {
00542       // delete old stuff if outa date
00543       if (theMaterials != 0) {
00544  for (int i=0; i<numFibers; i++)
00545    delete theMaterials[i];
00546  delete [] theMaterials;
00547  if (matData != 0)
00548    delete [] matData;
00549  matData = 0;
00550  theMaterials = 0;
00551       }
00552 
00553       // create memory to hold material pointers and fiber data
00554       numFibers = data(1);
00555       if (numFibers != 0) {
00556  theMaterials = new UniaxialMaterial *[numFibers];
00557  
00558  if (theMaterials == 0)
00559    g3ErrorHandler->fatal("%s -- failed to allocate Material pointers",
00560     "FiberSection2d::recvSelf");
00561  
00562  for (int j=0; j<numFibers; j++)
00563    theMaterials[j] = 0;
00564 
00565  matData = new double [numFibers*2];
00566 
00567  if (matData == 0)
00568    g3ErrorHandler->fatal("%s -- failed to allocate double array for material data",
00569     "FiberSection2d::recvSelf");
00570  
00571       }
00572     }
00573 
00574     Vector fiberData(matData, 3*numFibers);
00575     res += theChannel.recvVector(dbTag, commitTag, fiberData);
00576     if (res < 0) {
00577       g3ErrorHandler->warning("%s - failed to send material data",
00578          "FiberSection2d::recvSelf");
00579       return res;
00580     }    
00581 
00582     for (int i=0; i<numFibers; i++) {
00583       int classTag = materialData(2*i);
00584       int dbTag = materialData(2*i+1);
00585 
00586       // if material pointed to is blank or not of corrcet type, 
00587       // release old and create a new one
00588       if (theMaterials[i] == 0)
00589  theMaterials[i] = theBroker.getNewUniaxialMaterial(classTag);
00590       else if (theMaterials[i]->getClassTag() != classTag) {
00591  delete theMaterials[i];
00592  theMaterials[i] = theBroker.getNewUniaxialMaterial(classTag);      
00593       }
00594 
00595       if (theMaterials[i] == 0) 
00596  g3ErrorHandler->fatal("%s -- failed to allocate double array for material data",
00597          "FiberSection2d::recvSelf"); 
00598 
00599       theMaterials[i]->setDbTag(dbTag);
00600       res += theMaterials[i]->recvSelf(commitTag, theChannel, theBroker);
00601     }
00602   }    
00603 
00604 
00605   return res;
00606 }
00607 
00608 void
00609 FiberSection3d::Print(ostream &s, int flag)
00610 {
00611   s << "\nFiberSection3d, tag: " << this->getTag() << endl;
00612   s << "\tSection code: " << code;
00613   s << "\tNumber of Fibers: " << numFibers << endl;
00614   
00615   if (flag == 1)
00616    for (int i = 0; i < numFibers; i++)
00617     theMaterials[i]->Print(s, flag);
00618 }
00619 
00620 Response*
00621 FiberSection3d::setResponse(char **argv, int argc, Information &sectInfo)
00622 {
00623  // See if the response is one of the defaults
00624  Response *res = SectionForceDeformation::setResponse(argv, argc, sectInfo);
00625  if (res != 0)
00626   return res;
00627 
00628  // Check if fiber response is requested
00629     else if (strcmp(argv[0],"fiber") == 0) {
00630         int key = 0;
00631         int passarg = 2;
00632 
00633         if (argc <= 2)          // not enough data input
00634             return 0;
00635         else if (argc <= 3)  // fiber number was input directly
00636             key = atoi(argv[1]);
00637         else {                  // fiber near-to coordinate specified
00638             double yCoord = atof(argv[1]);
00639    double zCoord = atof(argv[2]);
00640    double ySearch = -matData[0];
00641    double zSearch =  matData[1];
00642    double closestDist = sqrt( pow(ySearch-yCoord,2) +
00643                                  pow(zSearch-zCoord,2) );
00644    double distance;
00645    for (int j = 1; j < numFibers; j++) {
00646     ySearch = -matData[3*j];
00647     zSearch =  matData[3*j+1];
00648     distance = sqrt( pow(ySearch-yCoord,2) +
00649                                  pow(zSearch-zCoord,2) );
00650     if (distance < closestDist) {
00651      closestDist = distance;
00652      key = j;
00653     }
00654    }
00655       passarg = 3;
00656   }
00657  
00658         if (key < numFibers)
00659    return theMaterials[key]->setResponse(&argv[passarg],argc-passarg,sectInfo);
00660         else
00661             return 0;
00662     }
00663    
00664     // otherwise response quantity is unknown for the FiberSection class
00665     else
00666   return 0;
00667 }
00668 
00669 
00670 int 
00671 FiberSection3d::getResponse(int responseID, Information &sectInfo)
00672 {
00673   // Just call the base class method ... don't need to define
00674   // this function, but keeping it here just for clarity
00675   return SectionForceDeformation::getResponse(responseID, sectInfo);
00676 }
00677