00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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
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
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
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
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
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
00114
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
00148
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
00181
00182 const Vector&
00183 FeapMaterial::getStress(void)
00184 {
00185 int isw = 3;
00186
00187
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
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
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
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
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
00376 for (i = 0; i < numHV; i++)
00377 vecData(i) = hstv[i];
00378
00379
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
00419 for (i = 0; i < numHV; i++)
00420 hstv[i] = vecData(i);
00421
00422
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
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
00501 double trace = eps[0] + eps[1] + eps[2];
00502
00503
00504 double td = 0.0;
00505
00506
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
00516 this->fillDArray();
00517
00518
00519 double dt = ops_Dt;
00520 int niter = 1;
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
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
00549 return 0;
00550 }