FeapMaterial.cpp

Go to the documentation of this file.
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.4 $
00022 // $Date: 2003/04/02 22:02:40 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/material/nD/FeapMaterial.cpp,v $
00024                                                                         
00025 // Written: MHS
00026 // Created: June 2001
00027 //
00028 // Description: This file contains the class definition for 
00029 // FeapMaterial. FeapMaterial wraps a Feap material subroutine.
00030 
00031 #include <OPS_Globals.h>
00032 #include <FeapMaterial.h>
00033 #include <Vector.h>
00034 #include <Matrix.h>
00035 #include <ID.h>
00036 #include <Channel.h>
00037 #include <string.h>
00038 #include <stdlib.h>
00039 #include <float.h>
00040 
00041 double FeapMaterial::d[200];
00042 double FeapMaterial::sig[6];
00043 double FeapMaterial::dd[36];
00044 
00045 Vector FeapMaterial::strain3(3);
00046 Vector FeapMaterial::strain4(4);
00047 Vector FeapMaterial::strain6(6);
00048 
00049 Vector FeapMaterial::sigma3(3);
00050 Vector FeapMaterial::sigma4(4);
00051 Vector FeapMaterial::sigma6(sig,6);
00052 
00053 Matrix FeapMaterial::tangent3(3,3);
00054 Matrix FeapMaterial::tangent4(4,4);
00055 Matrix FeapMaterial::tangent6(dd,6,6);
00056 
00057 FeapMaterial::FeapMaterial(int tag, int classTag, int nhv, int ndata, double r)
00058   :NDMaterial(tag,classTag), ud(0), hstv(0), rho(r),
00059    numHV(nhv), numData(ndata), myFormulation(ThreeDimensional)
00060 {
00061   if (numHV < 0)
00062     numHV = 0;
00063   
00064   if (numHV > 0) {
00065     // Allocate history array
00066     hstv = new double[2*numHV];
00067     if (hstv == 0) {
00068       opserr <<  "FeapMaterial::FeapMaterial -- failed to allocate history array -- type: " << classTag << endln;
00069       exit(-1);
00070     }
00071                             
00072     
00073     // Zero out the history variables
00074     for (int i = 0; i < 2*numHV; i++)
00075       hstv[i] = 0.0;
00076   }
00077   
00078   if (numData < 0)
00079     numData = 0;
00080   
00081   if (numData > 0) {
00082     // Allocate material parameter array
00083     ud = new double[numData];
00084     if (ud == 0) {
00085       opserr << "FeapMaterial::FeapMaterial -- failed to allocate ud array -- type: " << classTag << endln;
00086       exit(-1);
00087     }
00088   }
00089   
00090   // Zero the strain vector
00091   for (int i = 0; i < 6; i++)
00092     eps[i] = 0.0;
00093 }
00094 
00095 FeapMaterial::FeapMaterial(int classTag)
00096   :NDMaterial(0,classTag), ud(0), hstv(0), rho(0.0),
00097    numHV(0), numData(0), myFormulation(Unknown)
00098 {
00099   // Zero the strain vector
00100   for (int i = 0; i < 6; i++)
00101     eps[i] = 0.0;
00102 }
00103 
00104 FeapMaterial::~FeapMaterial()
00105 {
00106   if (ud != 0)
00107     delete [] ud;
00108   
00109   if (hstv != 0)
00110     delete [] hstv;
00111 }
00112 
00113 // Here we are assuming the following order on strains
00114 // \epsilon = {11, 22, 33, 12, 23, 31}
00115 int
00116 FeapMaterial::setTrialStrain(const Vector &strain)
00117 {
00118   switch (myFormulation) {
00119   case ThreeDimensional:
00120     eps[0] = strain(0);
00121     eps[1] = strain(1);
00122     eps[2] = strain(2);
00123     eps[3] = strain(3);
00124     eps[4] = strain(4);
00125     eps[5] = strain(5);
00126     break;
00127   case PlaneStrain:
00128     eps[0] = strain(0);
00129     eps[1] = strain(1);
00130     eps[3] = strain(2);
00131     break;
00132   case AxiSymmetric:
00133     eps[0] = strain(0);
00134     eps[1] = strain(1);
00135     eps[2] = strain(2);
00136     eps[3] = strain(3);
00137     break;
00138   default:
00139     opserr << "FeapMaterial::FeapMaterial -- unknown material formulation\n";
00140     exit(-1);
00141     break;
00142   }
00143 
00144   return 0;
00145 }
00146 
00147 // Here we are assuming the following order on strains
00148 // \epsilon = {11, 22, 33, 12, 23, 31}
00149 const Vector&
00150 FeapMaterial::getStrain(void)
00151 {
00152   switch (myFormulation) {
00153   case ThreeDimensional:
00154     strain6(0) = eps[0];
00155     strain6(1) = eps[1];
00156     strain6(2) = eps[2];
00157     strain6(3) = eps[3];
00158     strain6(4) = eps[4];
00159     strain6(5) = eps[5];
00160     return strain6;
00161   case PlaneStrain:
00162     strain3(0) = eps[0];
00163     strain3(1) = eps[1];
00164     strain3(2) = eps[3];
00165     return strain3;
00166   case AxiSymmetric:
00167     strain4(0) = eps[0];
00168     strain4(1) = eps[1];
00169     strain4(2) = eps[2];
00170     strain4(3) = eps[3];
00171     return strain4;
00172   default:
00173     opserr << "FeapMaterial::getSTrain -- unknown material formulation\n";
00174     exit(-1);
00175       
00176     return strain6;
00177   }
00178 }
00179 
00180 // Here we are assuming the following order on stresses
00181 // \sigma = {11, 22, 33, 12, 23, 31}
00182 const Vector&
00183 FeapMaterial::getStress(void)
00184 {
00185   int isw = 3;
00186 
00187   // Invoke Feap subroutine
00188   this->invokeSubroutine(isw);
00189   
00190   switch (myFormulation) {
00191   case ThreeDimensional:
00192     return sigma6;
00193   case PlaneStrain:
00194     sigma3(0) = sig[0];
00195     sigma3(1) = sig[1];
00196     sigma3(2) = sig[3];
00197     return sigma3;
00198   case AxiSymmetric:
00199     sigma4(0) = sig[0];
00200     sigma4(1) = sig[1];
00201     sigma4(2) = sig[2];
00202     sigma4(3) = sig[3];
00203     return sigma4;
00204   default:
00205     opserr << "FeapMaterial::getStress -- unknown material formulation\n";
00206       exit(-1);
00207     return sigma6;
00208   }
00209 }
00210 
00211 const Matrix&
00212 FeapMaterial::getTangent(void)
00213 {
00214   int isw = 6;
00215 
00216   // Invoke Feap subroutine
00217   this->invokeSubroutine(isw);
00218 
00219   int i,j;
00220 
00221   switch (myFormulation) {
00222   case ThreeDimensional:
00223     return tangent6;
00224   case PlaneStrain:
00225     tangent3(0,0) = tangent6(0,0);
00226     tangent3(0,1) = tangent6(0,1);
00227     tangent3(0,2) = tangent6(0,3);
00228     tangent3(1,0) = tangent6(1,0);
00229     tangent3(1,1) = tangent6(1,1);
00230     tangent3(1,2) = tangent6(1,3);
00231     tangent3(2,0) = tangent6(3,0);
00232     tangent3(2,1) = tangent6(3,1);
00233     tangent3(2,2) = tangent6(3,3);
00234     return tangent3;
00235   case AxiSymmetric:
00236     for (i = 0; i < 4; i++)
00237       for (j = 0; j < 4; j++)
00238         tangent4(i,j) = tangent6(i,j);
00239     return tangent4;
00240   default:
00241     opserr << "FeapMaterial::getTangent -- unknown material formulation\n";
00242     exit(-1);
00243     return tangent6;
00244   }
00245 }
00246 
00247 double
00248 FeapMaterial::getRho(void)
00249 {
00250   return rho;
00251 }
00252 
00253 int
00254 FeapMaterial::commitState(void)
00255 {
00256   // Set committed values equal to corresponding trial values
00257   for (int i = 0; i < numHV; i++)
00258     hstv[i] = hstv[i+numHV];
00259   
00260   return 0;
00261 }
00262 
00263 int
00264 FeapMaterial::revertToLastCommit(void)
00265 {
00266   // Set trial values equal to corresponding committed values
00267   for (int i = 0; i < numHV; i++)
00268     hstv[i+numHV] = hstv[i];
00269   
00270   return 0;
00271 }
00272 
00273 int
00274 FeapMaterial::revertToStart(void)
00275 {
00276   // Set all trial and committed values to zero
00277   for (int i = 0; i < 2*numHV; i++)
00278     hstv[i] = 0.0;
00279   
00280   return 0;
00281 }
00282 
00283 NDMaterial*
00284 FeapMaterial::getCopy(void)
00285 {
00286   FeapMaterial *theCopy = 
00287     new FeapMaterial(this->getTag(), this->getClassTag(),
00288                      numHV, numData, rho);
00289   
00290   int i;
00291   for (i = 0; i < 2*numHV; i++)
00292     theCopy->hstv[i] = hstv[i];
00293   
00294   for (i = 0; i < numData; i++)
00295     theCopy->ud[i] = ud[i];
00296   
00297   theCopy->myFormulation = myFormulation;
00298   
00299   return theCopy;
00300 }
00301 
00302 NDMaterial*
00303 FeapMaterial::getCopy(const char *type)
00304 {
00305   FeapMaterial *theCopy = (FeapMaterial*)this->getCopy();
00306   
00307   if (strcmp(type, "ThreeDimensional") == 0 || strcmp(type, "3D") == 0)
00308     theCopy->myFormulation = ThreeDimensional;
00309   
00310   else if (strcmp(type, "PlaneStrain") == 0 || strcmp(type, "PlaneStrain2D") == 0)
00311     theCopy->myFormulation = PlaneStrain;
00312   
00313   else if (strcmp(type, "AxiSymmetric") == 0 || strcmp(type, "AxiSymmetric2D") == 0)
00314     theCopy->myFormulation = AxiSymmetric;
00315   
00316   else {
00317     opserr << "FeapMaterial::getCopy -- Invalid type (" << type << ") for FeapMaterial\n";
00318     return 0;
00319   }
00320   
00321   return theCopy;
00322 }
00323 
00324 const char*
00325 FeapMaterial::getType(void) const
00326 {
00327   switch (myFormulation) {
00328   case ThreeDimensional:
00329     return "ThreeDimensional";
00330   case PlaneStrain:
00331     return "PlaneStrain";
00332   case AxiSymmetric:
00333     return "AxiSymmetric";
00334   default:
00335     opserr << "FeapMaterial::getTYpe -- unknown material formulation\n";
00336     return "Unknown";
00337   }
00338 }
00339 
00340 int
00341 FeapMaterial::getOrder(void) const
00342 {
00343   switch (myFormulation) {
00344   case ThreeDimensional:
00345     return 6;
00346   case PlaneStrain:
00347     return 3;
00348   case AxiSymmetric:
00349     return 4;
00350   default:
00351     opserr << "FeapMaterial::getOrder -- unknown material formulation\n";
00352     return 0;
00353   }
00354 }
00355 
00356 int 
00357 FeapMaterial::sendSelf(int commitTag, Channel &theChannel)
00358 {
00359   int res = 0;
00360   
00361   static ID idData(4);
00362   
00363   idData(0) = this->getTag();
00364   idData(1) = numHV;
00365   idData(2) = numData;
00366   idData(3) = myFormulation;
00367   
00368   res += theChannel.sendID(this->getDbTag(), commitTag, idData);
00369   if (res < 0) 
00370     opserr << "FeapMaterial::sendSelf() - failed to send ID data\n";
00371   
00372   Vector vecData(numHV+numData);
00373   
00374   int i, j;
00375   // Copy history variables into vector
00376   for (i = 0; i < numHV; i++)
00377     vecData(i) = hstv[i];
00378   
00379   // Copy material properties into vector
00380   for (i = 0, j = numHV; i < numData; i++, j++)
00381     vecData(j) = ud[i];
00382   
00383   res += theChannel.sendVector(this->getDbTag(), commitTag, vecData);
00384   if (res < 0) 
00385     opserr << "FeapMaterial::sendSelf() - failed to send Vector data\n";
00386   
00387   return res;
00388 }
00389 
00390 int
00391 FeapMaterial::recvSelf(int commitTag, Channel &theChannel,
00392                        FEM_ObjectBroker &theBroker)
00393 {
00394   int res = 0;
00395   
00396   static ID idData(4);
00397   
00398   res += theChannel.recvID(this->getDbTag(), commitTag, idData);
00399   if (res < 0) {
00400     opserr << "FeapMaterial::recvSelf() - failed to receive ID data\n";
00401     return res;
00402   }
00403   
00404   this->setTag(idData(0));
00405   numHV   = idData(1);
00406   numData = idData(2);
00407   myFormulation  = idData(3);
00408   
00409   Vector vecData(numHV+numData);
00410   
00411   res += theChannel.recvVector(this->getDbTag(), commitTag, vecData);
00412   if (res < 0) {
00413     opserr << "FeapMaterial::recvSelf() - failed to receive Vector data\n";
00414     return res;
00415   }
00416   
00417   int i, j;
00418   // Copy history variables from vector
00419   for (i = 0; i < numHV; i++)
00420     hstv[i] = vecData(i);
00421   
00422   // Copy material properties from vector
00423   for (i = 0, j = numHV; i < numData; i++, j++)
00424     ud[i] = vecData(j);
00425   
00426   return res;
00427 }
00428     
00429 void
00430 FeapMaterial::Print(OPS_Stream &s, int flag)
00431 {
00432   s << "FeapMaterial, tag: " << this->getTag() << endln;
00433   s << "Material formulation: " << this->getType() << endln;
00434   s << "Material subroutine: ";
00435   
00436   switch (this->getClassTag()) {
00437   case ND_TAG_FeapMaterial01:
00438     s << "matl01" << endln;
00439     break;
00440 
00441   case ND_TAG_FeapMaterial02:
00442     s << "matl02" << endln;
00443     break;
00444 
00445   case ND_TAG_FeapMaterial03:
00446     s << "matl03" << endln;
00447     break;
00448     
00449     // Add more cases as needed
00450   default:
00451     s << this->getClassTag() << endln;
00452     break;
00453   }
00454   
00455   s << "Material density: " << rho << endln;
00456 }
00457 
00458 #ifdef _WIN32
00459 
00460 extern "C" int FEAPCOMMON(double *dt, int *niter);
00461 
00462 extern "C" int MATL01(double *eps, double *trace, double *td, double *d,
00463                                double *ud, double *hn, double *h1, int *nh,
00464                                double *sig, double *dd, int *isw);
00465 
00466 extern "C" int MATL02(double *eps, double *trace, double *td, double *d,
00467                                double *ud, double *hn, double *h1, int *nh,
00468                                double *sig, double *dd, int *isw);
00469 
00470 extern "C" int MATL03(double *eps, double *trace, double *td, double *d,
00471                                double *ud, double *hn, double *h1, int *nh,
00472                                double *sig, double *dd, int *isw);
00473 
00474 #define feapcommon_ FEAPCOMMON
00475 #define matl01_ MATL01
00476 #define matl02_ MATL02
00477 #define matl03_ MATL03
00478 
00479 #else
00480 
00481 extern "C" int feapcommon_(double *dt, int *niter);
00482 
00483 extern "C" int matl01_(double *eps, double *trace, double *td, double *d,
00484                        double *ud, double *hn, double *h1, int *nh,
00485                        double *sig, double *dd, int *isw);
00486 
00487 extern "C" int matl02_(double *eps, double *trace, double *td, double *d,
00488                        double *ud, double *hn, double *h1, int *nh,
00489                        double *sig, double *dd, int *isw);
00490 
00491 extern "C" int matl03_(double *eps, double *trace, double *td, double *d,
00492                        double *ud, double *hn, double *h1, int *nh,
00493                        double *sig, double *dd, int *isw);
00494 
00495 #endif
00496 
00497 int
00498 FeapMaterial::invokeSubroutine(int isw)
00499 {
00500   // Trace of strain vector
00501   double trace = eps[0] + eps[1] + eps[2];
00502 
00503   // Temperature change (currently not used)
00504   double td = 0.0;
00505 
00506   // Zero out stress and tangent arrays as required to do so
00507   int i;
00508   for (i = 0; i < 6; i++) {
00509     sig[i] = 0.0;
00510     dd[i] = 0.0;
00511   }
00512   for ( ; i < 36; i++)
00513     dd[i] = 0.0;
00514 
00515   // Populate the FEAP material array for this model
00516   this->fillDArray();
00517 
00518   // Fill in the common blocks
00519   double dt = ops_Dt;  // From G3Globals.h
00520   int niter = 1;    // Need to count the number of global iterations!
00521   feapcommon_(&dt, &niter);
00522 
00523   switch (this->getClassTag()) {
00524   case ND_TAG_FeapMaterial01:
00525     matl01_(eps, &trace, &td, d, ud, hstv, &hstv[numHV], &numHV, sig, dd, &isw);
00526     break;
00527 
00528   case ND_TAG_FeapMaterial02:
00529     matl02_(eps, &trace, &td, d, ud, hstv, &hstv[numHV], &numHV, sig, dd, &isw);
00530     break;
00531 
00532   case ND_TAG_FeapMaterial03:
00533     matl03_(eps, &trace, &td, d, ud, hstv, &hstv[numHV], &numHV, sig, dd, &isw);
00534     break;
00535     
00536     // Add more cases as needed
00537   default:
00538     opserr << "FeapMaterial::invokeSubroutine -- unknown material type\n";
00539     return -1;
00540   }
00541   
00542   return 0;
00543 }
00544 
00545 int
00546 FeapMaterial::fillDArray(void)
00547 {
00548   // Must be done by subclasses!
00549   return 0;
00550 }

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