Template3Dep.cppGo to the documentation of this file.00001 /* 00002 00003 ################################################################################ 00004 00005 # COPYRIGHT (C): :-)) # 00006 00007 # PROJECT: Object Oriented Finite Element Program # 00008 00009 # PURPOSE: General platform for elaso-plastic constitutive model # 00010 00011 # implementation # 00012 00013 # CLASS: Template3Dep (the base class for all material point) # 00014 00015 # # 00016 00017 # VERSION: # 00018 00019 # LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) # 00020 00021 # TARGET OS: DOS || UNIX || . . . # 00022 00023 # DESIGNER(S): Boris Jeremic, Zhaohui Yang # 00024 00025 # PROGRAMMER(S): Boris Jeremic, Zhaohui Yang # 00026 00027 # # 00028 00029 # # 00030 00031 # DATE: 08-03-2000 # 00032 00033 # UPDATE HISTORY: 09-12-2000 # 00034 00035 # May 2004, Zhao Cheng splitting the elastic part # 00036 00037 # Oct. 2004 Zhao Cheng, small addition for u-p-U modeling # 00038 00039 # Mar. 2005 Guanzhou updated constitutive driver to be # 00040 00041 # compatible with global Newton-Raphson iterations# # 00042 00043 # BackwardEuler has been corrected # 00044 00045 # # 00046 00047 # SHORT EXPLANATION: This file contains the class implementation for # 00048 00049 # Template3Dep. # 00050 00051 # # 00052 00053 ################################################################################ 00054 00055 */ 00056 00057 00058 00059 #ifndef Template3Dep_CPP 00060 00061 #define Template3Dep_CPP 00062 00063 00064 00065 #define ITMAX 30 00066 00067 #define MAX_STEP_COUNT 30 00068 00069 #define NUM_OF_SUB_INCR 30 00070 00071 #define KK 1000.0 00072 00073 00074 00075 #include "Template3Dep.h" 00076 00077 00078 00079 00080 00081 stresstensor Template3Dep::Stress; 00082 00083 straintensor Template3Dep::Strain; 00084 00085 00086 00087 00088 00089 //================================================================================ 00090 00091 // Constructor 00092 00093 //================================================================================ 00094 00095 00096 00097 Template3Dep::Template3Dep( int tag , 00098 00099 NDMaterial &theElMat, 00100 00101 YieldSurface *YS_ , 00102 00103 PotentialSurface *PS_ , 00104 00105 EPState *EPS_ , 00106 00107 EvolutionLaw_S *ELS1_ , 00108 00109 EvolutionLaw_S *ELS2_ , 00110 00111 EvolutionLaw_S *ELS3_ , 00112 00113 EvolutionLaw_S *ELS4_ , 00114 00115 EvolutionLaw_T *ELT1_ , 00116 00117 EvolutionLaw_T *ELT2_ , 00118 00119 EvolutionLaw_T *ELT3_ , 00120 00121 EvolutionLaw_T *ELT4_ ) 00122 00123 :NDMaterial(tag, ND_TAG_Template3Dep) 00124 00125 { 00126 00127 // Get a copy of the material 00128 00129 theElasticMat = theElMat.getCopy(); 00130 00131 if (theElasticMat == 0) { 00132 00133 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00134 00135 exit(-1); 00136 00137 } 00138 00139 00140 00141 if ( YS_ ) 00142 00143 YS = YS_->newObj(); 00144 00145 else 00146 00147 { 00148 00149 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00150 00151 exit(-1); 00152 00153 } 00154 00155 00156 00157 if ( PS_ ) 00158 00159 PS = PS_->newObj(); 00160 00161 else 00162 00163 { 00164 00165 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00166 00167 exit(-1); 00168 00169 } 00170 00171 00172 00173 if ( EPS_ ) 00174 00175 EPS = EPS_->newObj(); 00176 00177 else 00178 00179 { 00180 00181 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00182 00183 exit(-1); 00184 00185 } 00186 00187 00188 00189 // Evolution laws 00190 00191 if ( ELS1_ ) 00192 00193 ELS1 = ELS1_->newObj(); 00194 00195 else 00196 00197 ELS1 = 0; 00198 00199 00200 00201 if ( ELS2_ ) 00202 00203 ELS2 = ELS2_->newObj(); 00204 00205 else 00206 00207 ELS2 = 0; 00208 00209 00210 00211 if ( ELS3_ ) 00212 00213 ELS3 = ELS3_->newObj(); 00214 00215 else 00216 00217 ELS3 = 0; 00218 00219 00220 00221 if ( ELS4_ ) 00222 00223 ELS4 = ELS4_->newObj(); 00224 00225 else 00226 00227 ELS4 = 0; 00228 00229 00230 00231 if ( ELT1_ ) 00232 00233 ELT1 = ELT1_->newObj(); 00234 00235 else 00236 00237 ELT1 = 0; 00238 00239 00240 00241 if ( ELT2_ ) 00242 00243 ELT2 = ELT2_->newObj(); 00244 00245 else 00246 00247 ELT2 = 0; 00248 00249 00250 00251 if ( ELT3_ ) 00252 00253 ELT3 = ELT3_->newObj(); 00254 00255 else 00256 00257 ELT3 = 0; 00258 00259 00260 00261 if ( ELT4_ ) 00262 00263 ELT4 = ELT4_->newObj(); 00264 00265 else 00266 00267 ELT4 = 0; 00268 00269 00270 00271 //Initialze Eep using E-elastic 00272 00273 tensor E = ElasticStiffnessTensor(); 00274 00275 EPS->setEep( E ); 00276 00277 00278 00279 } 00280 00281 00282 00283 //================================================================================ 00284 00285 // Constructor 0 00286 00287 //================================================================================ 00288 00289 Template3Dep::Template3Dep( int tag , 00290 00291 NDMaterial &theElMat, 00292 00293 YieldSurface *YS_ , 00294 00295 PotentialSurface *PS_ , 00296 00297 EPState *EPS_) 00298 00299 :NDMaterial(tag, ND_TAG_Template3Dep) 00300 00301 { 00302 00303 // Get a copy of the material 00304 00305 theElasticMat = theElMat.getCopy(); 00306 00307 if (theElasticMat == 0) { 00308 00309 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00310 00311 exit(-1); 00312 00313 } 00314 00315 00316 00317 if ( YS_ ) 00318 00319 YS = YS_->newObj(); 00320 00321 else 00322 00323 { 00324 00325 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00326 00327 exit(-1); 00328 00329 } 00330 00331 00332 00333 if ( PS_ ) 00334 00335 PS = PS_->newObj(); 00336 00337 else 00338 00339 { 00340 00341 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00342 00343 exit(-1); 00344 00345 } 00346 00347 00348 00349 if ( EPS_ ) 00350 00351 EPS = EPS_->newObj(); 00352 00353 else 00354 00355 { 00356 00357 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00358 00359 exit(-1); 00360 00361 } 00362 00363 00364 00365 // Evolution laws 00366 00367 ELS1 = 0; 00368 00369 ELS2 = 0; 00370 00371 ELS3 = 0; 00372 00373 ELS4 = 0; 00374 00375 ELT1 = 0; 00376 00377 ELT2 = 0; 00378 00379 ELT3 = 0; 00380 00381 ELT4 = 0; 00382 00383 } 00384 00385 00386 00387 00388 00389 //================================================================================ 00390 00391 // Constructor 1 00392 00393 //================================================================================ 00394 00395 00396 00397 Template3Dep::Template3Dep( int tag , 00398 00399 NDMaterial &theElMat, 00400 00401 YieldSurface *YS_ , 00402 00403 PotentialSurface *PS_ , 00404 00405 EPState *EPS_, 00406 00407 EvolutionLaw_S *ELS1_ ) 00408 00409 :NDMaterial(tag, ND_TAG_Template3Dep) 00410 00411 { 00412 00413 // Get a copy of the material 00414 00415 theElasticMat = theElMat.getCopy(); 00416 00417 if (theElasticMat == 0) { 00418 00419 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00420 00421 exit(-1); 00422 00423 } 00424 00425 00426 00427 if ( YS_ ) 00428 00429 YS = YS_->newObj(); 00430 00431 else 00432 00433 { 00434 00435 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00436 00437 exit(-1); 00438 00439 } 00440 00441 00442 00443 if ( PS_ ) 00444 00445 PS = PS_->newObj(); 00446 00447 else 00448 00449 { 00450 00451 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00452 00453 exit(-1); 00454 00455 } 00456 00457 00458 00459 if ( EPS_ ) 00460 00461 EPS = EPS_->newObj(); 00462 00463 else 00464 00465 { 00466 00467 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00468 00469 exit(-1); 00470 00471 } 00472 00473 00474 00475 // Evolution laws 00476 00477 if ( ELS1_ ) 00478 00479 ELS1 = ELS1_->newObj(); 00480 00481 else 00482 00483 ELS1 = 0; 00484 00485 00486 00487 ELS2 = 0; 00488 00489 ELS3 = 0; 00490 00491 ELS4 = 0; 00492 00493 ELT1 = 0; 00494 00495 ELT2 = 0; 00496 00497 ELT3 = 0; 00498 00499 ELT4 = 0; 00500 00501 } 00502 00503 00504 00505 00506 00507 //================================================================================ 00508 00509 // Constructor 2 00510 00511 //================================================================================ 00512 00513 00514 00515 Template3Dep::Template3Dep( int tag , 00516 00517 NDMaterial &theElMat, 00518 00519 YieldSurface *YS_ , 00520 00521 PotentialSurface *PS_ , 00522 00523 EPState *EPS_, 00524 00525 EvolutionLaw_T *ELT1_ ) 00526 00527 :NDMaterial(tag, ND_TAG_Template3Dep) 00528 00529 { 00530 00531 // Get a copy of the material 00532 00533 theElasticMat = theElMat.getCopy(); 00534 00535 if (theElasticMat == 0) { 00536 00537 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00538 00539 exit(-1); 00540 00541 } 00542 00543 00544 00545 if ( YS_ ) 00546 00547 YS = YS_->newObj(); 00548 00549 else 00550 00551 { 00552 00553 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00554 00555 exit(-1); 00556 00557 } 00558 00559 00560 00561 if ( PS_ ) 00562 00563 PS = PS_->newObj(); 00564 00565 else 00566 00567 { 00568 00569 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00570 00571 exit(-1); 00572 00573 } 00574 00575 00576 00577 if ( EPS_ ) 00578 00579 EPS = EPS_->newObj(); 00580 00581 else 00582 00583 { 00584 00585 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00586 00587 exit(-1); 00588 00589 } 00590 00591 00592 00593 // Evolution laws 00594 00595 ELS1 = 0; 00596 00597 ELS2 = 0; 00598 00599 ELS3 = 0; 00600 00601 ELS4 = 0; 00602 00603 00604 00605 if ( ELT1_ ) 00606 00607 ELT1 = ELT1_->newObj(); 00608 00609 else 00610 00611 ELT1 = 0; 00612 00613 00614 00615 ELT2 = 0; 00616 00617 ELT3 = 0; 00618 00619 ELT4 = 0; 00620 00621 } 00622 00623 00624 00625 //================================================================================ 00626 00627 // Constructor 3 00628 00629 //================================================================================ 00630 00631 00632 00633 Template3Dep::Template3Dep( int tag, 00634 00635 NDMaterial &theElMat, 00636 00637 YieldSurface *YS_ , 00638 00639 PotentialSurface *PS_ , 00640 00641 EPState *EPS_, 00642 00643 EvolutionLaw_S *ELS1_, 00644 00645 EvolutionLaw_T *ELT1_ ) 00646 00647 :NDMaterial(tag, ND_TAG_Template3Dep) 00648 00649 { 00650 00651 // Get a copy of the material 00652 00653 theElasticMat = theElMat.getCopy(); 00654 00655 if (theElasticMat == 0) { 00656 00657 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00658 00659 exit(-1); 00660 00661 } 00662 00663 00664 00665 if ( YS_ ) 00666 00667 YS = YS_->newObj(); 00668 00669 else 00670 00671 { 00672 00673 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00674 00675 exit(-1); 00676 00677 } 00678 00679 00680 00681 if ( PS_ ) 00682 00683 PS = PS_->newObj(); 00684 00685 else 00686 00687 { 00688 00689 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00690 00691 exit(-1); 00692 00693 } 00694 00695 00696 00697 if ( EPS_ ) 00698 00699 EPS = EPS_->newObj(); 00700 00701 else 00702 00703 { 00704 00705 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00706 00707 exit(-1); 00708 00709 } 00710 00711 00712 00713 // Evolution laws 00714 00715 if ( ELS1_ ) 00716 00717 ELS1 = ELS1_->newObj(); 00718 00719 else 00720 00721 ELS1 = 0; 00722 00723 ELS2 = 0; 00724 00725 ELS3 = 0; 00726 00727 ELS4 = 0; 00728 00729 00730 00731 if ( ELT1_ ) 00732 00733 ELT1 = ELT1_->newObj(); 00734 00735 else 00736 00737 ELT1 = 0; 00738 00739 ELT2 = 0; 00740 00741 ELT3 = 0; 00742 00743 ELT4 = 0; 00744 00745 } 00746 00747 00748 00749 //================================================================================ 00750 00751 // Constructor 4 00752 00753 // Two scalar evolution laws and one tensorial evolution law are provided! 00754 00755 //================================================================================ 00756 00757 00758 00759 Template3Dep::Template3Dep( int tag , 00760 00761 NDMaterial &theElMat, 00762 00763 YieldSurface *YS_ , 00764 00765 PotentialSurface *PS_ , 00766 00767 EPState *EPS_, 00768 00769 EvolutionLaw_S *ELS1_, 00770 00771 EvolutionLaw_S *ELS2_, 00772 00773 EvolutionLaw_T *ELT1_ ) 00774 00775 :NDMaterial(tag, ND_TAG_Template3Dep) 00776 00777 { 00778 00779 // Get a copy of the material 00780 00781 theElasticMat = theElMat.getCopy(); 00782 00783 if (theElasticMat == 0) { 00784 00785 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00786 00787 exit(-1); 00788 00789 } 00790 00791 00792 00793 if ( YS_ ) 00794 00795 YS = YS_->newObj(); 00796 00797 else 00798 00799 { 00800 00801 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00802 00803 exit(-1); 00804 00805 } 00806 00807 00808 00809 if ( PS_ ) 00810 00811 PS = PS_->newObj(); 00812 00813 else 00814 00815 { 00816 00817 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00818 00819 exit(-1); 00820 00821 } 00822 00823 00824 00825 if ( EPS_ ) 00826 00827 EPS = EPS_->newObj(); 00828 00829 else 00830 00831 { 00832 00833 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00834 00835 exit(-1); 00836 00837 } 00838 00839 00840 00841 if ( ELS1_ ) 00842 00843 ELS1 = ELS1_->newObj(); 00844 00845 else 00846 00847 ELS1 = 0; 00848 00849 00850 00851 if ( ELS2_ ) 00852 00853 ELS2 = ELS2_->newObj(); 00854 00855 else 00856 00857 ELS2 = 0; 00858 00859 00860 00861 ELS3 = 0; 00862 00863 ELS4 = 0; 00864 00865 00866 00867 if ( ELT1_ ) 00868 00869 ELT1 = ELT1_->newObj(); 00870 00871 else 00872 00873 ELT1 = 0; 00874 00875 00876 00877 ELT2 = 0; 00878 00879 ELT3 = 0; 00880 00881 ELT4 = 0; 00882 00883 } 00884 00885 00886 00887 //================================================================================ 00888 00889 // Constructor 5 00890 00891 // Two scalar evolution laws and two tensorial evolution law are provided! 00892 00893 //================================================================================ 00894 00895 00896 00897 Template3Dep::Template3Dep( int tag , 00898 00899 NDMaterial &theElMat, 00900 00901 YieldSurface *YS_ , 00902 00903 PotentialSurface *PS_ , 00904 00905 EPState *EPS_, 00906 00907 EvolutionLaw_S *ELS1_, 00908 00909 EvolutionLaw_S *ELS2_, 00910 00911 EvolutionLaw_T *ELT1_, 00912 00913 EvolutionLaw_T *ELT2_) 00914 00915 :NDMaterial(tag, ND_TAG_Template3Dep) 00916 00917 { 00918 00919 // Get a copy of the material 00920 00921 theElasticMat = theElMat.getCopy(); 00922 00923 if (theElasticMat == 0) { 00924 00925 opserr << "Template3Dep:: Template3Dep failed to get copy of material\n"; 00926 00927 exit(-1); 00928 00929 } 00930 00931 00932 00933 if ( YS_ ) 00934 00935 YS = YS_->newObj(); 00936 00937 else 00938 00939 { 00940 00941 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00942 00943 exit(-1); 00944 00945 } 00946 00947 00948 00949 if ( PS_ ) 00950 00951 PS = PS_->newObj(); 00952 00953 else 00954 00955 { 00956 00957 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00958 00959 exit(-1); 00960 00961 } 00962 00963 00964 00965 if ( EPS_ ) 00966 00967 EPS = EPS_->newObj(); 00968 00969 else 00970 00971 { 00972 00973 opserr << "Template3Dep:: Template3Dep failed to construct the template3Dep\n"; 00974 00975 exit(-1); 00976 00977 } 00978 00979 00980 00981 if ( ELS1_ ) 00982 00983 ELS1 = ELS1_->newObj(); 00984 00985 else 00986 00987 ELS1 = 0; 00988 00989 00990 00991 if ( ELS2_ ) 00992 00993 ELS2 = ELS2_->newObj(); 00994 00995 else 00996 00997 ELS2 = 0; 00998 00999 ELS3 = 0; 01000 01001 ELS4 = 0; 01002 01003 01004 01005 if ( ELT1_ ) 01006 01007 ELT1 = ELT1_->newObj(); 01008 01009 else 01010 01011 ELT1 = 0; 01012 01013 01014 01015 if ( ELT2_ ) 01016 01017 ELT2 = ELT2_->newObj(); 01018 01019 else 01020 01021 ELT2 = 0; 01022 01023 ELT3 = 0; 01024 01025 ELT4 = 0; 01026 01027 } 01028 01029 01030 01031 //================================================================================ 01032 01033 // Constructor 6 01034 01035 // NO parameter is provided 01036 01037 //================================================================================ 01038 01039 01040 01041 Template3Dep::Template3Dep() 01042 01043 :NDMaterial(0, ND_TAG_Template3Dep),ELS1(0), ELS2(0),ELS3(0), ELS4(0), 01044 01045 ELT1(0), ELT2(0), ELT3(0), ELT4(0) 01046 01047 { 01048 01049 //YS = new DPYieldSurface(); 01050 01051 //PS = new DPPotentialSurface(); 01052 01053 //EPS = new EPState(); 01054 01055 theElasticMat = 0; 01056 01057 YS = 0; 01058 01059 PS = 0; 01060 01061 EPS = 0; 01062 01063 } 01064 01065 01066 01067 01068 01069 //================================================================================ 01070 01071 // Destructor 01072 01073 //================================================================================ 01074 01075 01076 01077 Template3Dep::~Template3Dep() 01078 01079 { 01080 01081 if (theElasticMat) 01082 01083 delete theElasticMat; 01084 01085 01086 01087 if (YS) 01088 01089 delete YS; 01090 01091 01092 01093 if (PS) 01094 01095 delete PS; 01096 01097 01098 01099 if (EPS) 01100 01101 delete EPS; 01102 01103 01104 01105 if (ELS1) 01106 01107 delete ELS1; 01108 01109 01110 01111 if (ELS2) 01112 01113 delete ELS2; 01114 01115 01116 01117 if (ELS3) 01118 01119 delete ELS3; 01120 01121 01122 01123 if (ELS4) 01124 01125 delete ELS4; 01126 01127 01128 01129 if (ELT1) 01130 01131 delete ELT1; 01132 01133 01134 01135 if (ELT2) 01136 01137 delete ELT2; 01138 01139 01140 01141 if (ELT3) 01142 01143 delete ELT3; 01144 01145 01146 01147 if (ELT4) 01148 01149 delete ELT4; 01150 01151 01152 01153 01154 01155 } 01156 01157 01158 01160 01162 01164 01165 //Template3Dep::Template3Dep(const Template3Dep & rval) { 01166 01167 // 01168 01169 // YS = rval.YS->newObj(); 01170 01171 // PS = rval.PS->newObj(); 01172 01173 // EPS = rval.EPS->newObj(); 01174 01175 // 01176 01177 // // Scalar Evolution Laws 01178 01179 // if ( rval.getELS1() ) 01180 01181 // ELS1 = rval.getELS1()->newObj(); 01182 01183 // else 01184 01185 // ELS1 = 0; 01186 01187 // 01188 01189 // if ( !rval.getELS2() ) 01190 01191 // ELS2 = 0; 01192 01193 // else 01194 01195 // ELS2 = rval.getELS2()->newObj(); 01196 01197 // 01198 01199 // if ( !rval.getELS3() ) 01200 01201 // ELS3 = 0; 01202 01203 // else 01204 01205 // ELS3 = rval.getELS3()->newObj(); 01206 01207 // 01208 01209 // if ( !rval.getELS4() ) 01210 01211 // ELS4 = 0; 01212 01213 // else 01214 01215 // ELS4 = rval.getELS4()->newObj(); 01216 01217 // 01218 01219 // // Tensorial Evolution Laws 01220 01221 // if ( rval.getELT1() ) 01222 01223 // ELT1 = rval.getELT1()->newObj(); 01224 01225 // else 01226 01227 // ELT1 = 0; 01228 01229 // 01230 01231 // if ( !rval.getELT2() ) 01232 01233 // ELT2 = 0; 01234 01235 // else 01236 01237 // ELT2 = rval.getELT2()->newObj(); 01238 01239 // 01240 01241 // if ( !rval.getELT3() ) 01242 01243 // ELT3 = 0; 01244 01245 // else 01246 01247 // ELT3 = rval.getELT3()->newObj(); 01248 01249 // 01250 01251 // if ( !rval.getELT4() ) 01252 01253 // ELT4 = 0; 01254 01255 // else 01256 01257 // ELT4 = rval.getELT4()->newObj(); 01258 01259 // 01260 01261 //} 01262 01263 // 01264 01265 01266 01267 01268 01269 01270 01271 //================================================================================ 01272 01273 // Routine used to generate elastic compliance tensor D for this material point 01274 01275 //================================================================================ 01276 01277 tensor Template3Dep::ElasticComplianceTensor(void) const 01278 01279 { 01280 01281 tensor ret( 4, def_dim_4, 0.0 ); 01282 01283 01284 01285 ret = ElasticStiffnessTensor().inverse(); 01286 01287 01288 01289 return ret; 01290 01291 01292 01293 //ZC05/2004 int eflag = (this->EPS)->getElasticflag(); 01294 01295 //ZC05/2004 double Ey = this->EPS->getEo(); 01296 01297 //ZC05/2004 double nuP =this->EPS->getnu(); 01298 01299 //ZC05/2004 01300 01301 //ZC05/2004 // Zhao, 03-09-04 01302 01303 //ZC05/2004 // eflag = 0, Pressure independent Isotropic 01304 01305 //ZC05/2004 // = 1, Pressure dendent Isotropic 01306 01307 //ZC05/2004 // = 2, Pressure independent Cross Anisotropic 01308 01309 //ZC05/2004 // = 3, Pressure dependent Cross Anisotropic 01310 01311 //ZC05/2004 // default: 0, Pressure indendent Isotropic 01312 01313 //ZC05/2004 01314 01315 //ZC05/2004 // Pressure independent Isotropic 01316 01317 //ZC05/2004 if ( eflag == 0 ) 01318 01319 //ZC05/2004 { 01320 01321 //ZC05/2004 if (Ey == 0) 01322 01323 //ZC05/2004 { 01324 01325 //ZC05/2004 opserr << endln << "Ey = 0! Can't give you D!!" << endln; 01326 01327 //ZC05/2004 exit(1); 01328 01329 //ZC05/2004 } 01330 01331 //ZC05/2004 01332 01333 //ZC05/2004 // Kronecker delta tensor 01334 01335 //ZC05/2004 tensor I2("I", 2, def_dim_2); 01336 01337 //ZC05/2004 01338 01339 //ZC05/2004 tensor I_ijkl = I2("ij")*I2("kl"); 01340 01341 //ZC05/2004 //I_ijkl.null_indices(); 01342 01343 //ZC05/2004 tensor I_ikjl = I_ijkl.transpose0110(); 01344 01345 //ZC05/2004 tensor I_iljk = I_ijkl.transpose0111(); 01346 01347 //ZC05/2004 tensor I4s = (I_ikjl+I_iljk)*0.5; 01348 01349 //ZC05/2004 01350 01351 //ZC05/2004 // Building compliance tensor 01352 01353 //ZC05/2004 ret = (I_ijkl * (-nuP/Ey)) + (I4s * ((1.0+nuP)/Ey)); 01354 01355 //ZC05/2004 } 01356 01357 //ZC05/2004 01358 01359 //ZC05/2004 // Pressure dependent Isotropic 01360 01361 //ZC05/2004 else if ( eflag == 1) 01362 01363 //ZC05/2004 { 01364 01365 //ZC05/2004 double E = Ey; 01366 01367 //ZC05/2004 //opserr << " Eo= " << Ey; 01368 01369 //ZC05/2004 if (Ey == 0) 01370 01371 //ZC05/2004 { 01372 01373 //ZC05/2004 opserr << endln << "Ey = 0! Can't give you D!!" << endln; 01374 01375 //ZC05/2004 exit(1); 01376 01377 //ZC05/2004 } 01378 01379 //ZC05/2004 01380 01381 //ZC05/2004 //Update E 01382 01383 //ZC05/2004 stresstensor stc = (this->EPS)->getStress(); 01384 01385 //ZC05/2004 double p = stc.p_hydrostatic(); 01386 01387 //ZC05/2004 double po = EPS->getpo(); 01388 01389 //ZC05/2004 01390 01391 //ZC05/2004 //opserr << " p = " << p; 01392 01393 //ZC05/2004 01394 01395 //ZC05/2004 //double po = 100.0; //kPa 01396 01397 //ZC05/2004 //BJ: this is questionable to be re-examined!!!!!!!!!!!!!!!!!!!! 01398 01399 //ZC05/2004 if (p <= 0.5000*KK) p = 0.500*KK; 01400 01401 //ZC05/2004 E = Ey * pow(p/po/KK, 0.5); //0.5 01402 01403 //ZC05/2004 //opserr << " Ec = " << Ey << endlnn; 01404 01405 //ZC05/2004 01406 01407 //ZC05/2004 // Kronecker delta tensor 01408 01409 //ZC05/2004 tensor I2("I", 2, def_dim_2); 01410 01411 //ZC05/2004 01412 01413 //ZC05/2004 tensor I_ijkl = I2("ij")*I2("kl"); 01414 01415 //ZC05/2004 //I_ijkl.null_indices(); 01416 01417 //ZC05/2004 tensor I_ikjl = I_ijkl.transpose0110(); 01418 01419 //ZC05/2004 tensor I_iljk = I_ijkl.transpose0111(); 01420 01421 //ZC05/2004 tensor I4s = (I_ikjl+I_iljk)*0.5; 01422 01423 //ZC05/2004 01424 01425 //ZC05/2004 // Building compliance tensor 01426 01427 //ZC05/2004 ret = (I_ijkl * (-nuP/E)) + (I4s * ((1.0+nuP)/E)); 01428 01429 //ZC05/2004 } 01430 01431 //ZC05/2004 01432 01433 //ZC05/2004 // Pressure independent Anisotropic 01434 01435 //ZC05/2004 else if (eflag == 2) 01436 01437 //ZC05/2004 { 01438 01439 //ZC05/2004 // Form compliance matrix D 01440 01441 //ZC05/2004 double nu = nuP; 01442 01443 //ZC05/2004 double Ev = this->EPS->getEv(); 01444 01445 //ZC05/2004 double Ghv = this->EPS->getGhv(); 01446 01447 //ZC05/2004 double nuhv = this->EPS->getnuhv(); 01448 01449 //ZC05/2004 01450 01451 //ZC05/2004 double A = 1.0/Ey; 01452 01453 //ZC05/2004 double B = 1.0/Ev; 01454 01455 //ZC05/2004 //opserr << "A " << A << " B " << B; 01456 01457 //ZC05/2004 01458 01459 //ZC05/2004 Matrix D(6,6); 01460 01461 //ZC05/2004 D(0,0) = D(1,1) = A; 01462 01463 //ZC05/2004 D(2,2) = B; 01464 01465 //ZC05/2004 D(0,1) = D(1,0) = -nu*A; 01466 01467 //ZC05/2004 D(0,2) = D(2,0) = D(1,2) = D(2,1) = -nuhv*B; 01468 01469 //ZC05/2004 //D(3,3) = (1.0+nu)*A; // Bug found, 03-09-04 01470 01471 //ZC05/2004 //D(4,4) = D(5,5) = 0.5/Ghv; // ... 01472 01473 //ZC05/2004 D(3,3) = 2.0*(1.0+nu)*A; 01474 01475 //ZC05/2004 D(4,4) = D(5,5) = 1.0/Ghv; 01476 01477 //ZC05/2004 //opserr << " C " << D; 01478 01479 //ZC05/2004 01480 01481 //ZC05/2004 //Convert Matric D to 4th order Elastic constants tensor ret; 01482 01483 //ZC05/2004 ret.val(1,1,1,1) = D(0,0); ret.val(1,1,2,2) = D(0,1); ret.val(1,1,3,3) = D(0,2); // --> Sigma_xx 01484 01485 //ZC05/2004 ret.val(1,2,1,2) = D(3,3); ret.val(1,2,2,1) = D(3,3); // --> Sigma_xy 01486 01487 //ZC05/2004 ret.val(1,3,1,3) = D(4,4); ret.val(1,3,3,1) = D(4,4); // --> Sigma_xz 01488 01489 //ZC05/2004 01490 01491 //ZC05/2004 ret.val(2,1,1,2) = D(3,3); ret.val(2,1,2,1) = D(3,3); // --> Sigma_yx 01492 01493 //ZC05/2004 ret.val(2,2,1,1) = D(1,0); ret.val(2,2,2,2) = D(1,1); ret.val(2,2,3,3) = D(1,2); // --> Sigma_yy 01494 01495 //ZC05/2004 ret.val(2,3,2,3) = D(5,5); ret.val(2,3,3,2) = D(5,5); // --> Sigma_yz 01496 01497 //ZC05/2004 01498 01499 //ZC05/2004 ret.val(3,1,1,3) = D(4,4); ret.val(3,1,3,1) = D(4,4); // --> Sigma_zx 01500 01501 //ZC05/2004 ret.val(3,2,2,3) = D(5,5); ret.val(3,2,3,2) = D(5,5); // --> Sigma_zy 01502 01503 //ZC05/2004 ret.val(3,3,1,1) = D(2,0); ret.val(3,3,2,2) = D(2,1); ret.val(3,3,3,3) = D(2,2); // --> Sigma_zz 01504 01505 //ZC05/2004 01506 01507 //ZC05/2004 } 01508 01509 //ZC05/2004 01510 01511 //ZC05/2004 // Pressure dependent Anisotropic 01512 01513 //ZC05/2004 else if (eflag == 3) 01514 01515 //ZC05/2004 { 01516 01517 //ZC05/2004 // Form compliance matrix D 01518 01519 //ZC05/2004 double E = Ey; 01520 01521 //ZC05/2004 double nu = nuP; 01522 01523 //ZC05/2004 double Ev = this->EPS->getEv(); 01524 01525 //ZC05/2004 double Ghv = this->EPS->getGhv(); 01526 01527 //ZC05/2004 double nuhv = this->EPS->getnuhv(); 01528 01529 //ZC05/2004 //opserr << " Eo= " << Ey; 01530 01531 //ZC05/2004 if (Ey == 0) 01532 01533 //ZC05/2004 { 01534 01535 //ZC05/2004 opserr << endln << "Ey = 0! Can't give you D!!" << endln; 01536 01537 //ZC05/2004 exit(1); 01538 01539 //ZC05/2004 } 01540 01541 //ZC05/2004 01542 01543 //ZC05/2004 //Update E 01544 01545 //ZC05/2004 stresstensor stc = (this->EPS)->getStress(); 01546 01547 //ZC05/2004 double p = stc.p_hydrostatic(); 01548 01549 //ZC05/2004 double po = EPS->getpo(); 01550 01551 //ZC05/2004 01552 01553 //ZC05/2004 //opserr << " p = " << p; 01554 01555 //ZC05/2004 01556 01557 //ZC05/2004 //double po = 100.0; //kPa 01558 01559 //ZC05/2004 if (p <= 0.5000*KK) p = 0.500*KK; 01560 01561 //ZC05/2004 Ey = Ey * pow(p/po/KK, 0.5); //0.5 01562 01563 //ZC05/2004 Ev = Ev * pow(p/po/KK, 0.5); //0.5 01564 01565 //ZC05/2004 Ghv = Ghv * pow(p/po/KK, 0.5); //0.5 01566 01567 //ZC05/2004 //opserr << " Ec = " << Ey << endlnn; 01568 01569 //ZC05/2004 01570 01571 //ZC05/2004 double A = 1.0/E; 01572 01573 //ZC05/2004 double B = 1.0/Ev; 01574 01575 //ZC05/2004 //opserr << "A " << A << " B " << B; 01576 01577 //ZC05/2004 01578 01579 //ZC05/2004 Matrix D(6,6); 01580 01581 //ZC05/2004 D(0,0) = D(1,1) = A; 01582 01583 //ZC05/2004 D(2,2) = B; 01584 01585 //ZC05/2004 D(0,1) = D(1,0) = -nu*A; 01586 01587 //ZC05/2004 D(0,2) = D(2,0) = D(1,2) = D(2,1) = -nuhv*B; 01588 01589 //ZC05/2004 D(3,3) = 2.0*(1.0+nu)*A; 01590 01591 //ZC05/2004 D(4,4) = D(5,5) = 1.0/Ghv; 01592 01593 //ZC05/2004 //opserr << " C " << D; 01594 01595 //ZC05/2004 01596 01597 //ZC05/2004 //Convert Matric D to 4th order Elastic constants tensor ret; 01598 01599 //ZC05/2004 ret.val(1,1,1,1) = D(0,0); ret.val(1,1,2,2) = D(0,1); ret.val(1,1,3,3) = D(0,2); // --> Sigma_xx 01600 01601 //ZC05/2004 ret.val(1,2,1,2) = D(3,3); ret.val(1,2,2,1) = D(3,3); // --> Sigma_xy 01602 01603 //ZC05/2004 ret.val(1,3,1,3) = D(4,4); ret.val(1,3,3,1) = D(4,4); // --> Sigma_xz 01604 01605 //ZC05/2004 01606 01607 //ZC05/2004 ret.val(2,1,1,2) = D(3,3); ret.val(2,1,2,1) = D(3,3); // --> Sigma_yx 01608 01609 //ZC05/2004 ret.val(2,2,1,1) = D(1,0); ret.val(2,2,2,2) = D(1,1); ret.val(2,2,3,3) = D(1,2); // --> Sigma_yy 01610 01611 //ZC05/2004 ret.val(2,3,2,3) = D(5,5); ret.val(2,3,3,2) = D(5,5); // --> Sigma_yz 01612 01613 //ZC05/2004 01614 01615 //ZC05/2004 ret.val(3,1,1,3) = D(4,4); ret.val(3,1,3,1) = D(4,4); // --> Sigma_zx 01616 01617 //ZC05/2004 ret.val(3,2,2,3) = D(5,5); ret.val(3,2,3,2) = D(5,5); // --> Sigma_zy 01618 01619 //ZC05/2004 ret.val(3,3,1,1) = D(2,0); ret.val(3,3,2,2) = D(2,1); ret.val(3,3,3,3) = D(2,2); // --> Sigma_zz 01620 01621 //ZC05/2004 01622 01623 //ZC05/2004 } 01624 01625 //ZC05/2004 else 01626 01627 //ZC05/2004 { 01628 01629 //ZC05/2004 opserr << "Template3Dep::ElasticComplianceTensor failed to create elastic compliance tensor. Eflag must be 0-3: " << eflag << endln; 01630 01631 //ZC05/2004 exit(-1); 01632 01633 //ZC05/2004 01634 01635 //ZC05/2004 } 01636 01637 //ZC05/2004 01638 01639 //ZC05/2004 return ret; 01640 01641 } 01642 01643 01644 01645 01646 01647 //================================================================================ 01648 01649 // Routine used to generate elastic stiffness tensor D for this material point 01650 01651 //================================================================================ 01652 01653 tensor Template3Dep::ElasticStiffnessTensor(void) const 01654 01655 { 01656 01657 tensor ret( 4, def_dim_4, 0.0 ); 01658 01659 01660 01661 if (theElasticMat->setTrialStrain( this->EPS->getElasticStrain() ) < 0) { 01662 01663 opserr << "Template3Dep::theElasticMat setTrialStrain failed in material with strain "; 01664 01665 return -1; 01666 01667 } 01668 01669 01670 01671 ret = theElasticMat->getTangentTensor(); 01672 01673 return ret; 01674 01675 01676 01677 //ZC05/2004 int eflag = (this->EPS)->getElasticflag(); 01678 01679 //ZC05/2004 01680 01681 //ZC05/2004 double Ey = this->EPS->getEo(); 01682 01683 //ZC05/2004 double nu =this->EPS->getnu(); 01684 01685 //ZC05/2004 01686 01687 //ZC05/2004 // Zhao, 03-09-04 01688 01689 //ZC05/2004 // eflag = 0, Pressure dependent Isotropic 01690 01691 //ZC05/2004 // = 1, Pressure indendent Isotropic 01692 01693 //ZC05/2004 // = 2, Pressure dependent Cross Anisotropic 01694 01695 //ZC05/2004 // = 3, Pressure independent Cross Anisotropic 01696 01697 //ZC05/2004 // default: 0, Pressure indendent Isotropic 01698 01699 //ZC05/2004 // Pressure dependent Isotropic 01700 01701 //ZC05/2004 if ( eflag == 0) { 01702 01703 //ZC05/2004 double E = Ey; 01704 01705 //ZC05/2004 // Kronecker delta tensor 01706 01707 //ZC05/2004 tensor I2("I", 2, def_dim_2); 01708 01709 //ZC05/2004 01710 01711 //ZC05/2004 tensor I_ijkl = I2("ij")*I2("kl"); 01712 01713 //ZC05/2004 01714 01715 //ZC05/2004 01716 01717 //ZC05/2004 //I_ijkl.null_indices(); 01718 01719 //ZC05/2004 tensor I_ikjl = I_ijkl.transpose0110(); 01720 01721 //ZC05/2004 tensor I_iljk = I_ijkl.transpose0111(); 01722 01723 //ZC05/2004 tensor I4s = (I_ikjl+I_iljk)*0.5; 01724 01725 //ZC05/2004 01726 01727 //ZC05/2004 // Building elasticity tensor 01728 01729 //ZC05/2004 ret = I_ijkl*( E*nu / ( (1.0+nu)*(1.0 - 2.0*nu) ) ) + I4s*( E / (1.0 + nu) ); 01730 01731 //ZC05/2004 //ret.null_indices(); 01732 01733 //ZC05/2004 } 01734 01735 //ZC05/2004 else if ( eflag == 1) { 01736 01737 //ZC05/2004 01738 01739 //ZC05/2004 //Update E 01740 01741 //ZC05/2004 stresstensor stc = (this->EPS)->getStress(); 01742 01743 //ZC05/2004 double p = stc.p_hydrostatic(); 01744 01745 //ZC05/2004 double po = EPS->getpo(); 01746 01747 //ZC05/2004 //opserr << " p = " << p; 01748 01749 //ZC05/2004 01750 01751 //ZC05/2004 //double po = 100.0; //kPa 01752 01753 //ZC05/2004 if (p <= 0.5000*KK) p = 0.5000*KK; 01754 01755 //ZC05/2004 double E = Ey * pow(p/po/KK, 0.5);//0.5 01756 01757 //ZC05/2004 //opserr << " Eo = " << Ey ; 01758 01759 //ZC05/2004 //opserr << " Ec = " << E << endlnn; 01760 01761 //ZC05/2004 01762 01763 //ZC05/2004 01764 01765 //ZC05/2004 // Kronecker delta tensor 01766 01767 //ZC05/2004 tensor I2("I", 2, def_dim_2); 01768 01769 //ZC05/2004 01770 01771 //ZC05/2004 tensor I_ijkl = I2("ij")*I2("kl"); 01772 01773 //ZC05/2004 01774 01775 //ZC05/2004 01776 01777 //ZC05/2004 //I_ijkl.null_indices(); 01778 01779 //ZC05/2004 tensor I_ikjl = I_ijkl.transpose0110(); 01780 01781 //ZC05/2004 tensor I_iljk = I_ijkl.transpose0111(); 01782 01783 //ZC05/2004 tensor I4s = (I_ikjl+I_iljk)*0.5; 01784 01785 //ZC05/2004 01786 01787 //ZC05/2004 // Building elasticity tensor 01788 01789 //ZC05/2004 ret = I_ijkl*( E*nu / ( (1.0+nu)*(1.0 - 2.0*nu) ) ) + I4s*( E / (1.0 + nu) ); 01790 01791 //ZC05/2004 //ret.null_indices(); 01792 01793 //ZC05/2004 } 01794 01795 //ZC05/2004 else if ( eflag == 3) { 01796 01797 //ZC05/2004 double E = Ey; 01798 01799 //ZC05/2004 double Ev = this->EPS->getEv(); 01800 01801 //ZC05/2004 double Ghv = this->EPS->getGhv(); 01802 01803 //ZC05/2004 double nuhv = this->EPS->getnuhv(); 01804 01805 //ZC05/2004 //Update E 01806 01807 //ZC05/2004 01808 01809 //ZC05/2004 stresstensor stc = (this->EPS)->getStress(); 01810 01811 //ZC05/2004 double p = stc.p_hydrostatic(); 01812 01813 //ZC05/2004 double po = EPS->getpo(); 01814 01815 //ZC05/2004 01816 01817 //ZC05/2004 //double po = 100.0; //kPa 01818 01819 //ZC05/2004 if (p <= 0.5000*KK) p = 0.500*KK; 01820 01821 //ZC05/2004 E = Ey * pow(p/po/KK, 0.5); //0.5 01822 01823 //ZC05/2004 Ev = Ev * pow(p/po/KK, 0.5); //0.5 01824 01825 //ZC05/2004 Ghv = Ghv * pow(p/po/KK, 0.5); //0.5 01826 01827 //ZC05/2004 //opserr << " Ec = " << Ey << endlnn; 01828 01829 //ZC05/2004 double A = 1.0/E; 01830 01831 //ZC05/2004 double B = 1.0/Ev; 01832 01833 //ZC05/2004 //opserr << "A " << A << " B " << B; 01834 01835 //ZC05/2004 01836 01837 //ZC05/2004 Matrix D(6,6); 01838 01839 //ZC05/2004 D(0,0) = D(1,1) = A; 01840 01841 //ZC05/2004 D(2,2) = B; 01842 01843 //ZC05/2004 D(0,1) = D(1,0) = -nu*A; 01844 01845 //ZC05/2004 D(0,2) = D(2,0) = D(1,2) = D(2,1) = -nuhv*B; 01846 01847 //ZC05/2004 D(3,3) = (1.0+nu)*A; 01848 01849 //ZC05/2004 D(4,4) = D(5,5) = 0.5/Ghv; 01850 01851 //ZC05/2004 //opserr << " C " << D; 01852 01853 //ZC05/2004 01854 01855 //ZC05/2004 // Do invertion once to get Elastic matrix D 01856 01857 //ZC05/2004 D.Invert( D ); 01858 01859 //ZC05/2004 01860 01861 //ZC05/2004 //Convert Matric D to 4th order Elastic constants tensor ret; 01862 01863 //ZC05/2004 ret.val(1,1,1,1) = D(0,0); ret.val(1,1,2,2) = D(0,1); ret.val(1,1,3,3) = D(0,2); // --> Sigma_xx 01864 01865 //ZC05/2004 ret.val(1,2,1,2) = D(3,3); ret.val(1,2,2,1) = D(3,3); // --> Sigma_xy 01866 01867 //ZC05/2004 ret.val(1,3,1,3) = D(4,4); ret.val(1,3,3,1) = D(4,4); // --> Sigma_xz 01868 01869 //ZC05/2004 01870 01871 //ZC05/2004 ret.val(2,1,1,2) = D(3,3); ret.val(2,1,2,1) = D(3,3); // --> Sigma_yx 01872 01873 //ZC05/2004 ret.val(2,2,1,1) = D(1,0); ret.val(2,2,2,2) = D(1,1); ret.val(2,2,3,3) = D(1,2); // --> Sigma_yy 01874 01875 //ZC05/2004 ret.val(2,3,2,3) = D(5,5); ret.val(2,3,3,2) = D(5,5); // --> Sigma_yz 01876 01877 //ZC05/2004 01878 01879 //ZC05/2004 ret.val(3,1,1,3) = D(4,4); ret.val(3,1,3,1) = D(4,4); // --> Sigma_zx 01880 01881 //ZC05/2004 ret.val(3,2,2,3) = D(5,5); ret.val(3,2,3,2) = D(5,5); // --> Sigma_zy 01882 01883 //ZC05/2004 ret.val(3,3,1,1) = D(2,0); ret.val(3,3,2,2) = D(2,1); ret.val(3,3,3,3) = D(2,2); // --> Sigma_zz 01884 01885 //ZC05/2004 01886 01887 //ZC05/2004 } 01888 01889 //ZC05/2004 else if ( eflag == 2) { 01890 01891 //ZC05/2004 double Ev = this->EPS->getEv(); 01892 01893 //ZC05/2004 double Ghv = this->EPS->getGhv(); 01894 01895 //ZC05/2004 double nuhv = this->EPS->getnuhv(); 01896 01897 //ZC05/2004 01898 01899 //ZC05/2004 double A = 1.0/Ey; 01900 01901 //ZC05/2004 double B = 1.0/Ev; 01902 01903 //ZC05/2004 //opserr << "A " << A << " B " << B; 01904 01905 //ZC05/2004 01906 01907 //ZC05/2004 Matrix D(6,6); 01908 01909 //ZC05/2004 D(0,0) = D(1,1) = A; 01910 01911 //ZC05/2004 D(2,2) = B; 01912 01913 //ZC05/2004 D(0,1) = D(1,0) = -nu*A; 01914 01915 //ZC05/2004 D(0,2) = D(2,0) = D(1,2) = D(2,1) = -nuhv*B; 01916 01917 //ZC05/2004 D(3,3) = (1.0+nu)*A; 01918 01919 //ZC05/2004 D(4,4) = D(5,5) = 0.5/Ghv; 01920 01921 //ZC05/2004 //opserr << " C " << D; 01922 01923 //ZC05/2004 01924 01925 //ZC05/2004 // Do invertion once to get Elastic matrix D 01926 01927 //ZC05/2004 D.Invert( D ); 01928 01929 //ZC05/2004 01930 01931 //ZC05/2004 //Convert Matric D to 4th order Elastic constants tensor ret; 01932 01933 //ZC05/2004 ret.val(1,1,1,1) = D(0,0); ret.val(1,1,2,2) = D(0,1); ret.val(1,1,3,3) = D(0,2); // --> Sigma_xx 01934 01935 //ZC05/2004 ret.val(1,2,1,2) = D(3,3); ret.val(1,2,2,1) = D(3,3); // --> Sigma_xy 01936 01937 //ZC05/2004 ret.val(1,3,1,3) = D(4,4); ret.val(1,3,3,1) = D(4,4); // --> Sigma_xz 01938 01939 //ZC05/2004 01940 01941 //ZC05/2004 ret.val(2,1,1,2) = D(3,3); ret.val(2,1,2,1) = D(3,3); // --> Sigma_yx 01942 01943 //ZC05/2004 ret.val(2,2,1,1) = D(1,0); ret.val(2,2,2,2) = D(1,1); ret.val(2,2,3,3) = D(1,2); // --> Sigma_yy 01944 01945 //ZC05/2004 ret.val(2,3,2,3) = D(5,5); ret.val(2,3,3,2) = D(5,5); // --> Sigma_yz 01946 01947 //ZC05/2004 01948 01949 //ZC05/2004 ret.val(3,1,1,3) = D(4,4); ret.val(3,1,3,1) = D(4,4); // --> Sigma_zx 01950 01951 //ZC05/2004 ret.val(3,2,2,3) = D(5,5); ret.val(3,2,3,2) = D(5,5); // --> Sigma_zy 01952 01953 //ZC05/2004 ret.val(3,3,1,1) = D(2,0); ret.val(3,3,2,2) = D(2,1); ret.val(3,3,3,3) = D(2,2); // --> Sigma_zz 01954 01955 //ZC05/2004 01956 01957 //ZC05/2004 } 01958 01959 //ZC05/2004 else 01960 01961 //ZC05/2004 { 01962 01963 //ZC05/2004 opserr << "Template3Dep::ElasticStiffnessTensor failed to create elastic compliance tensor. Eflag must be 0-3: " << eflag << endln; 01964 01965 //ZC05/2004 exit(-1); 01966 01967 //ZC05/2004 01968 01969 //ZC05/2004 } 01970 01971 //ZC05/2004 return ret; 01972 01973 } 01974 01975 01976 01977 //================================================================================ 01978 01979 int Template3Dep::setTrialStrain(const Vector &v) 01980 01981 { 01982 01983 // Not yet implemented 01984 01985 return 0; 01986 01987 } 01988 01989 01990 01991 //================================================================================ 01992 01993 int Template3Dep::setTrialStrain(const Vector &v, const Vector &r) 01994 01995 { 01996 01997 // Not yet implemented 01998 01999 return 0; 02000 02001 } 02002 02003 02004 02005 //================================================================================ 02006 02007 int Template3Dep::setTrialStrainIncr(const Vector &v) 02008 02009 { 02010 02011 // Not yet implemented 02012 02013 return 0; 02014 02015 } 02016 02017 02018 02019 //================================================================================ 02020 02021 int Template3Dep::setTrialStrainIncr(const Vector &v, const Vector &r) 02022 02023 { 02024 02025 //================================================================================ 02026 02027 // Not yet implemented 02028 02029 return 0; 02030 02031 } 02032 02033 02034 02035 //================================================================================ 02036 02037 const Matrix& Template3Dep::getTangent(void) 02038 02039 { 02040 02041 // Not yet implemented 02042 02043 Matrix *M = new Matrix(); 02044 02045 return *M; 02046 02047 } 02048 02049 02050 02051 const Matrix& Template3Dep::getInitialTangent(void) 02052 02053 { 02054 02055 return this->getInitialTangent(); 02056 02057 } 02058 02059 02060 02061 //================================================================================ 02062 02063 const Vector& Template3Dep::getStress(void) 02064 02065 { 02066 02067 // Not yet implemented 02068 02069 Vector *V = new Vector(); 02070 02071 return *V; 02072 02073 } 02074 02075 02076 02077 //================================================================================ 02078 02079 const Vector& Template3Dep::getStrain(void) 02080 02081 { 02082 02083 // Not yet implemented 02084 02085 Vector *V = new Vector(); 02086 02087 return *V; 02088 02089 } 02090 02091 02092 02093 //================================================================================ 02094 02095 // what is the trial strain? Initial strain? 02096 02097 int Template3Dep::setTrialStrain(const Tensor &v) 02098 02099 { 02100 02101 //Guanzhou made it compatible with global iterations Mar2005 02102 02103 this->setTrialStrainIncr( v - EPS->getStrain() ); 02104 02105 return 0; 02106 02107 } 02108 02109 02110 02111 02112 02113 //================================================================================ 02114 02115 int Template3Dep::setTrialStrain(const Tensor &v, const Tensor &r) 02116 02117 { 02118 02119 //Guanzhou made it compatible with global iterations Mar2005 02120 02121 this->setTrialStrainIncr( v - EPS->getStrain() ); 02122 02123 return 0; 02124 02125 } 02126 02127 02128 02129 //================================================================================ 02130 02131 02132 02133 int Template3Dep::setTrialStrainIncr(const Tensor &v) 02134 02135 { 02136 02137 02138 02139 //opserr << "\nBE: " << endlnn; 02140 02141 //EPState StartEPS( *(this->getEPS()) ); 02142 02143 //stresstensor start_stress = StartEPS.getStress_commit(); //Guanzhou Mar2005 02144 02145 //opserr << "start_stress 0 " << start_stress; 02146 02147 02148 02149 //EPState tmp_EPS = BackwardEulerEPState(v); 02150 02151 //if ( tmp_EPS.getConverged() ) { 02152 02153 // //setTrialEPS( tmp_EPS ); 02154 02155 // setEPS( tmp_EPS ); 02156 02157 // //double p = (tmp_EPS.getStress()).p_hydrostatic(); 02158 02159 // //double ec = (tmp_EPS.getec()) - (tmp_EPS.getLam()) * log( p / (tmp_EPS.getpo()) ); 02160 02161 // //double st = (tmp_EPS.getStrain()).Iinvariant1(); 02162 02163 // //double pl_s = (tmp_EPS.getPlasticStrain()).Iinvariant1(); 02164 02165 // //double dpl_s = (tmp_EPS.getdPlasticStrain()).Iinvariant1(); 02166 02167 // //opserr << "ec " << ec << " e " << tmp_EPS.gete() << " psi " << (tmp_EPS.gete() - ec) << " strain " << st << " t_pl " << pl_s << " d_pl " << dpl_s << "\n"; 02168 02169 // return 0; 02170 02171 //} 02172 02173 02174 02175 //opserr << endlnn; 02176 02177 //setEPS( StartEPS ); 02178 02179 //int number_of_subincrements = 5; 02180 02181 //Cascading subdividing in case that some incr_step is too big 02182 02183 02184 02185 //int loop = 0; 02186 02187 //while ( !tmp_EPS.getConverged() && (loop < 1) ) { 02188 02189 // 02190 02191 // setEPS( StartEPS ); 02192 02193 // EPState startEPS( *(this->getEPS()) ); 02194 02195 // stresstensor start_stress = startEPS.getStress(); 02196 02197 // opserr << " Step Start Stress:" << start_stress << endln; 02198 02199 // 02200 02201 // loop += 1; 02202 02203 // opserr << "\n "<< loop << "th Sub-BE, total subdivision: " << 10*loop*NUM_OF_SUB_INCR << endln; 02204 02205 // tmp_EPS = BESubIncrementation(v, 10*loop*NUM_OF_SUB_INCR); 02206 02207 // if ( tmp_EPS.getConverged() ) { 02208 02209 // //setTrialEPS( tmp_EPS ); 02210 02211 // setEPS( tmp_EPS ); 02212 02213 // return 0; 02214 02215 // } 02216 02217 // //else { 02218 02219 // // opserr << "\n2nd Sub BE: " << 3*NUM_OF_SUB_INCR << endln; 02220 02221 // // tmp_EPS = BESubIncrementation(v, 3*NUM_OF_SUB_INCR); 02222 02223 // // 02224 02225 // // if ( tmp_EPS.getConverged() ) { 02226 02227 // // setEPS( tmp_EPS ); 02228 02229 // // return 0; 02230 02231 // // } 02232 02233 // // else 02234 02235 // // return 1; 02236 02237 // //} 02238 02239 //} 02240 02241 02242 02244 02245 //EPState tmp_EPS = FESubIncrementation(v, NUM_OF_SUB_INCR); 02246 02247 EPState *thisEPState = this->getEPS(); 02248 02249 EPState tmp_EPS; 02250 02251 if ( thisEPState->getIntegratorFlag() == 0 ) tmp_EPS = ForwardEulerEPState(v); 02252 02253 else if ( thisEPState->getIntegratorFlag() == 1 ) tmp_EPS = BackwardEulerEPState(v); 02254 02255 else { 02256 02257 opserr << "Template3Dep::setTrialStrainIncr, error when getting integrator flag from EPState! \n"; 02258 02259 exit(1); 02260 02261 } 02262 02263 02264 02265 //EPState tmp_EPS = BackwardEulerEPState(v); 02266 02267 setEPS( tmp_EPS ); 02268 02269 //setEPS( StartEPS ); 02270 02271 return 0; 02272 02273 } 02274 02275 02276 02277 //================================================================================ 02278 02279 int Template3Dep::setTrialStrainIncr(const Tensor &v, const Tensor &r) 02280 02281 { 02282 02283 //EPS->setStrain(v + EPS->getStrain() ); //ZC10/26/2004 02284 02285 this->setTrialStrainIncr(v); //ZC10/26/2004 02286 02287 return 0; 02288 02289 } 02290 02291 02292 02293 //================================================================================ 02294 02295 const tensor& Template3Dep::getTangentTensor(void) 02296 02297 { 02298 02299 // Tensor Eep = EPS->getEep(); 02300 02301 // return Eep; 02302 02303 return EPS->Eep; 02304 02305 } 02306 02307 02308 02309 //================================================================================ 02310 02311 const stresstensor& Template3Dep::getStressTensor(void) 02312 02313 { 02314 02315 //opserr << *EPS; 02316 02317 //stresstensor tmp; 02318 02319 //tmp = EPS->getStress(); 02320 02321 //opserr << "EPS->getStress() " << EPS->getStress() << endlnn; 02322 02323 02324 02325 //Something funny!!! happened here when returning EPS->getStress() 02326 02327 // This function will return wrong Stress. 02328 02329 //stresstensor ret = EPS->getStress(); 02330 02331 //return ret; 02332 02333 Stress = EPS->getStress(); 02334 02335 return Stress; 02336 02337 } 02338 02339 02340 02341 02342 02343 //================================================================================ 02344 02345 const straintensor& Template3Dep::getStrainTensor(void) 02346 02347 { 02348 02349 Strain = EPS->getStrain(); 02350 02351 return Strain; 02352 02353 } 02354 02355 02356 02357 //================================================================================ 02358 02359 const straintensor& Template3Dep::getPlasticStrainTensor(void) 02360 02361 { 02362 02363 Strain = EPS->getPlasticStrain(); 02364 02365 return Strain; 02366 02367 } 02368 02369 02370 02372 02373 //double Template3Dep::getpsi(void) 02374 02375 //{ 02376 02377 // //this function cannot be moved, 02378 02379 // //leave here for compiling..., Zhao04/22/04 02380 02381 // return 0.05; 02382 02383 //} 02384 02385 02386 02387 //================================================================================ 02388 02389 int Template3Dep::commitState(void) 02390 02391 { 02392 02393 int err; 02394 02395 theElasticMat->commitState(); 02396 02397 err = getEPS()->commitState(); 02398 02399 return err; 02400 02401 } 02402 02403 02404 02405 //================================================================================ 02406 02407 int Template3Dep::revertToLastCommit(void) 02408 02409 { 02410 02411 int err; 02412 02413 theElasticMat->revertToLastCommit(); 02414 02415 err = EPS->revertToLastCommit(); 02416 02417 return err; 02418 02419 } 02420 02421 //================================================================================ 02422 02423 int Template3Dep::revertToStart(void) 02424 02425 { 02426 02427 int err; 02428 02429 theElasticMat->revertToStart(); 02430 02431 err = EPS->revertToStart(); 02432 02433 return err; 02434 02435 } 02436 02437 02438 02440 02442 02444 02445 //{ 02446 02447 // Template3Dep *t3d; 02448 02449 // return t3d; 02450 02451 //} 02452 02453 02454 02455 //================================================================================ 02456 02457 // Get one copy of the NDMaterial model 02458 02459 NDMaterial * Template3Dep::getCopy(void) 02460 02461 { 02462 02463 NDMaterial * tmp = 02464 02465 new Template3Dep( this->getTag() , 02466 02467 *(this->getElMat()), 02468 02469 this->getYS() , 02470 02471 this->getPS() , 02472 02473 this->getEPS() , 02474 02475 this->getELS1() , 02476 02477 this->getELS2() , 02478 02479 this->getELS3() , 02480 02481 this->getELS4() , 02482 02483 this->getELT1() , 02484 02485 this->getELT2() , 02486 02487 this->getELT3() , 02488 02489 this->getELT4() ); 02490 02491 02492 02493 return tmp; 02494 02495 02496 02497 } 02498 02499 02500 02501 02502 02503 //================================================================================ 02504 02505 NDMaterial * Template3Dep::getCopy(const char *code) 02506 02507 { 02508 02509 if (strcmp(code,"Template3Dep") == 0) 02510 02511 { 02512 02513 Template3Dep * tmp = 02514 02515 new Template3Dep( this->getTag() , 02516 02517 *(this->getElMat()), 02518 02519 this->getYS() , 02520 02521 this->getPS() , 02522 02523 this->getEPS() , 02524 02525 this->getELS1() , 02526 02527 this->getELS2() , 02528 02529 this->getELS3() , 02530 02531 this->getELS4() , 02532 02533 this->getELT1() , 02534 02535 this->getELT2() , 02536 02537 this->getELT3() , 02538 02539 this->getELT4() ); 02540 02541 02542 02543 tmp->getEPS()->setInit(); 02544 02545 return tmp; 02546 02547 } 02548 02549 else 02550 02551 { 02552 02553 opserr << "Template3Dep::getCopy failed to get model: " << code << endln; 02554 02555 exit(-1); 02556 02557 } 02558 02559 return 0; 02560 02561 } 02562 02563 02564 02565 //================================================================================ 02566 02567 const char *Template3Dep::getType(void) const 02568 02569 { 02570 02571 return "Template3Dep"; 02572 02573 } 02574 02575 02576 02577 //================================================================================ 02578 02579 //??What is the Order????????? might be the 02580 02581 // 02582 02583 //int Template3Dep::getOrder(void) const 02584 02585 //{ 02586 02587 // return 6; 02588 02589 //} 02590 02591 02592 02593 //================================================================================ 02594 02595 int Template3Dep::sendSelf(int commitTag, Channel &theChannel) 02596 02597 { 02598 02599 // Not yet implemented 02600 02601 return 0; 02602 02603 } 02604 02605 02606 02607 //================================================================================ 02608 02609 int Template3Dep::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) 02610 02611 { 02612 02613 // Not yet implemented 02614 02615 return 0; 02616 02617 } 02618 02619 02620 02621 //================================================================================ 02622 02623 void 02624 02625 Template3Dep::Print(OPS_Stream &s, int flag) 02626 02627 { 02628 02629 s << (*this); 02630 02631 } 02632 02633 02634 02635 02636 02637 02638 02639 //private utilities 02640 02641 02642 02643 //================================================================================ 02644 02645 // Set new EPState 02646 02647 //================================================================================ 02648 02649 02650 02651 void Template3Dep::setEPS( EPState & rval) 02652 02653 { 02654 02655 //EPState eps = rval; older buggy one 02656 02657 //EPS = rval.newObj(); 02658 02659 /* 02660 02661 //EPS->setEo(rval.getEo()); 02662 02663 EPS->setE(rval.getE()); 02664 02665 //EPS->setnu(getnu()); 02666 02667 //EPS->setrho(getrho()); 02668 02669 EPS->setStress(rval.getStress()); 02670 02671 EPS->setStrain(rval.getStrain()); 02672 02673 EPS->setElasticStrain(rval.getElasticStrain()); 02674 02675 EPS->setPlasticStrain(rval.getPlasticStrain()); 02676 02677 EPS->setdElasticStrain(rval.getdElasticStrain()); 02678 02679 EPS->setdPlasticStrain(rval.getdPlasticStrain()); 02680 02681 //EPS->setNScalarVar(rval.getNScalarVar()); 02682 02683 for (int i = 0; i <rval.getNScalarVar(); i++) 02684 02685 EPS->setScalarVar(i, rval.getScalarVar(i)); 02686 02687 02688 02689 02690 02691 //EPS->setNTensorVar(rval.getNTensorVar()); 02692 02693 EPS->setTensorVar(rval.getTensorVar()); 02694 02695 EPS->setEep(rval.getEep()); 02696 02697 EPS->Stress_commit=(rval.getStress_commit()); 02698 02699 EPS->Strain_commit=(rval.getStrain_commit()); 02700 02701 02702 02703 for (int i = 0; i <rval.getNScalarVar(); i++) 02704 02705 EPS->ScalarVar_commit[i] = rval.getScalarVar_commit()[i]; 02706 02707 02708 02709 for (int i = 0; i <rval.getNTensorVar(); i++) 02710 02711 EPS->TensorVar_commit[i] = (rval.getTensorVar_commit()[i]); 02712 02713 02714 02715 EPS->Eep_commit = (rval.getEep_commit()); 02716 02717 //EPS->setStress_init(getStress_init()); 02718 02719 //EPS->setStrain_init(getStrain_init()); 02720 02721 //EPS->setScalarVar_init(getScalarVar_init()); 02722 02723 //EPS->setTensorVar_init(getTensorVar_init()); 02724 02725 //EPS->setEep_init(getEep_init()); 02726 02727 EPS->setConverged(rval.getConverged()); 02728 02729 02730 02731 */ 02732 02733 02734 02735 //ZC05/2004 EPS->setElasticflag(rval.getElasticflag()); 02736 02737 //ZC05/2004 02738 02739 //ZC05/2004 EPS->setEo(rval.getEo()); 02740 02741 //ZC05/2004 EPS->setE(rval.getE()); 02742 02743 //ZC05/2004 EPS->setEv(rval.getEv()); 02744 02745 //ZC05/2004 EPS->setGhv(rval.getGhv()); 02746 02747 //ZC05/2004 EPS->setnu(rval.getnu()); 02748 02749 //ZC05/2004 EPS->setnuhv(rval.getnuhv()); 02750 02751 02752 02753 EPS->setStress(rval.getStress()); 02754 02755 EPS->setStrain(rval.getStrain()); 02756 02757 EPS->setElasticStrain(rval.getElasticStrain()); 02758 02759 EPS->setPlasticStrain(rval.getPlasticStrain()); 02760 02761 EPS->setdElasticStrain(rval.getdElasticStrain()); 02762 02763 EPS->setdPlasticStrain(rval.getdPlasticStrain()); 02764 02765 EPS->setNScalarVar( rval.getNScalarVar() ); 02766 02767 02768 02769 int i; 02770 02771 for (i = 0; i <rval.getNScalarVar(); i++) 02772 02773 EPS->setScalarVar(i+1, rval.getScalarVar(i+1)); 02774 02775 02776 02777 02778 02779 EPS->setNTensorVar( rval.getNTensorVar() ); 02780 02781 EPS->setTensorVar(rval.getTensorVar()); 02782 02783 EPS->setEep(rval.getEep()); 02784 02785 EPS->setStress_commit(rval.getStress_commit()); 02786 02787 EPS->setStrain_commit(rval.getStrain_commit()); 02788 02789 EPS->setElasticStrain_commit(rval.getElasticStrain_commit()); 02790 02791 02792 02793 for (i = 0; i <rval.getNScalarVar(); i++) 02794 02795 EPS->setScalarVar_commit(i+1, rval.getScalarVar_commit(i+1)); 02796 02797 02798 02799 for (i = 0; i <rval.getNTensorVar(); i++) 02800 02801 EPS->setTensorVar_commit(i+1, rval.getTensorVar_commit(i+1)); 02802 02803 02804 02805 EPS->Eep_commit = (rval.getEep_commit()); 02806 02807 EPS->Stress_init = rval.getStress_init(); 02808 02809 EPS->Strain_init = rval.getStrain_init(); 02810 02811 02812 02813 for (i = 0; i <rval.getNScalarVar(); i++) 02814 02815 EPS->setScalarVar_init(i+1, rval.getScalarVar_init(i+1)); 02816 02817 02818 02819 for (i = 0; i <rval.getNTensorVar(); i++) 02820 02821 EPS->setTensorVar_init(i+1, rval.getTensorVar_init(i+1)); 02822 02823 02824 02825 EPS->Eep_init = rval.getEep_init(); 02826 02827 EPS->setConverged(rval.getConverged()); 02828 02829 02830 02831 //ZC05/2004 //Added Joey 02-13-03 02832 02833 //ZC05/2004 EPS->seteo(rval.geteo()); 02834 02835 //ZC05/2004 EPS->setec(rval.getec()); 02836 02837 //ZC05/2004 EPS->setLam(rval.getLam()); 02838 02839 //ZC05/2004 EPS->setpo(rval.getpo()); 02840 02841 EPS->sete(rval.gete()); 02842 02843 EPS->setpsi(rval.getpsi()); 02844 02845 //ZC05/2004 EPS->seta(rval.geta()); 02846 02847 } 02848 02849 02850 02851 //================================================================================ 02852 02853 // Get the Elastic Material //ZC05/2004 02854 02855 //================================================================================ 02856 02857 NDMaterial* Template3Dep::getElMat() const 02858 02859 { 02860 02861 return theElasticMat; 02862 02863 } 02864 02865 02866 02867 02868 02869 //================================================================================ 02870 02871 // Get the Yield Surface 02872 02873 //================================================================================ 02874 02875 YieldSurface * Template3Dep::getYS() const 02876 02877 { 02878 02879 return YS; 02880 02881 } 02882 02883 02884 02885 02886 02887 //================================================================================ 02888 02889 // Get the Potential Surface 02890 02891 //================================================================================ 02892 02893 PotentialSurface * Template3Dep::getPS() const 02894 02895 { 02896 02897 return PS; 02898 02899 } 02900 02901 02902 02903 //================================================================================ 02904 02905 // Get the EPState 02906 02907 //================================================================================ 02908 02909 EPState * Template3Dep::getEPS() const 02910 02911 { 02912 02913 return EPS; 02914 02915 } 02916 02917 02918 02919 02920 02921 //================================================================================ 02922 02923 // Get the 1st Salar evolution law 02924 02925 //================================================================================ 02926 02927 02928 02929 EvolutionLaw_S * Template3Dep::getELS1() const 02930 02931 { 02932 02933 return ELS1; 02934 02935 } 02936 02937 02938 02939 //================================================================================ 02940 02941 // Get the 2ndst Salar evolution law 02942 02943 //================================================================================ 02944 02945 EvolutionLaw_S * Template3Dep::getELS2() const 02946 02947 { 02948 02949 return ELS2; 02950 02951 } 02952 02953 02954 02955 //================================================================================ 02956 02957 // Get the 2ndst Salar evolution law 02958 02959 //================================================================================ 02960 02961 EvolutionLaw_S * Template3Dep::getELS3() const 02962 02963 { 02964 02965 return ELS3; 02966 02967 } 02968 02969 //================================================================================ 02970 02971 // Get the 2ndst Salar evolution law 02972 02973 //================================================================================ 02974 02975 EvolutionLaw_S * Template3Dep::getELS4() const 02976 02977 { 02978 02979 return ELS4; 02980 02981 } 02982 02983 02984 02985 02986 02987 //================================================================================ 02988 02989 // Get the 1st tensorial evolution law 02990 02991 //================================================================================ 02992 02993 02994 02995 EvolutionLaw_T * Template3Dep::getELT1() const 02996 02997 { 02998 02999 return ELT1; 03000 03001 } 03002 03003 03004 03005 //================================================================================ 03006 03007 // Get the 2nd tensorial evolution law 03008 03009 //================================================================================ 03010 03011 03012 03013 EvolutionLaw_T * Template3Dep::getELT2() const 03014 03015 { 03016 03017 return ELT2; 03018 03019 } 03020 03021 //================================================================================ 03022 03023 // Get the 3rd tensorial evolution law 03024 03025 //================================================================================ 03026 03027 03028 03029 EvolutionLaw_T * Template3Dep::getELT3() const 03030 03031 { 03032 03033 return ELT3; 03034 03035 } 03036 03037 //================================================================================ 03038 03039 // Get the 4th tensorial evolution law 03040 03041 //================================================================================ 03042 03043 03044 03045 EvolutionLaw_T * Template3Dep::getELT4() const 03046 03047 { 03048 03049 return ELT4; 03050 03051 } 03052 03053 03054 03055 03056 03057 //================================================================================ 03058 03059 // Predictor EPState by Forward, Backward, MidPoint Methods... 03060 03061 //================================================================================ 03062 03063 03064 03065 EPState Template3Dep::PredictorEPState(straintensor & strain_increment) 03066 03067 { 03068 03069 03070 03071 EPState Predictor = ForwardEulerEPState( strain_increment); 03072 03073 //EPState Predictor = SemiBackwardEulerEPState( strain_increment, material_point); 03074 03075 return Predictor; 03076 03077 03078 03079 } 03080 03081 03082 03083 //================================================================================ 03084 03085 // New EPState using Forward Euler Algorithm 03086 03087 //================================================================================ 03088 03089 EPState Template3Dep::ForwardEulerEPState( const straintensor &strain_increment) 03090 03091 { 03092 03093 // Volumetric strain 03094 03095 //double st_vol = strain_increment.p_hydrostatic(); 03096 03097 double st_vol = strain_increment.Iinvariant1(); //- = compressive 03098 03099 //opserr << st_vol << " st_vol1 " << st_vol1 << "\n"; 03100 03101 03102 03103 //EPState forwardEPS( *(material_point->getEPS()) ); 03104 03105 EPState forwardEPS( *(this->getEPS()) ); 03106 03107 //opserr <<"start eps: " << forwardEPS; 03108 03109 //opserr << "\nForwardEulerEPState strain_increment " << strain_increment << endlnn; 03110 03111 03112 03113 // Building elasticity tensor 03114 03115 tensor E = ElasticStiffnessTensor(); 03116 03117 //tensor Eep = ElasticStiffnessTensor(); 03118 03119 tensor Eep = E; 03120 03121 tensor D = ElasticComplianceTensor(); 03122 03123 E.null_indices(); 03124 03125 D.null_indices(); 03126 03127 03128 03129 //Checking E and D 03130 03131 //tensor ed = E("ijpq") * D("pqkl"); 03132 03133 //double edd = ed.trace(); // = 3. 03134 03135 03136 03137 straintensor strain_incr = strain_increment; 03138 03139 strain_incr.null_indices(); 03140 03141 stresstensor stress_increment = E("ijpq") * strain_incr("pq"); 03142 03143 stress_increment.null_indices(); 03144 03145 //opserr << " stress_increment: " << stress_increment << endlnn; 03146 03147 03148 03149 EPState startEPS( *(getEPS()) ); 03150 03151 stresstensor start_stress = startEPS.getStress(); 03152 03153 start_stress.null_indices(); 03154 03155 //opserr << "===== start_EPS =====: " << startEPS; 03156 03157 03158 03159 double f_start = 0.0; 03160 03161 double f_pred = 0.0; 03162 03163 03164 03165 EPState IntersectionEPS( startEPS ); 03166 03167 03168 03169 EPState ElasticPredictorEPS( startEPS ); 03170 03171 stresstensor elastic_predictor_stress = start_stress + stress_increment; 03172 03173 ElasticPredictorEPS.setStress( elastic_predictor_stress ); 03174 03175 //opserr << " Elastic_predictor_stress: " << elastic_predictor_stress << endlnn; 03176 03177 03178 03179 f_start = this->getYS()->f( &startEPS ); 03180 03181 //::printf("\n############## f_start = %.10e ",f_start); 03182 03183 //opserr << "\nFE f_start = " << f_start; 03184 03185 03186 03187 f_pred = this->getYS()->f( &ElasticPredictorEPS ); 03188 03189 //::printf("############## f_pred = %.10e\n\n",f_pred); 03190 03191 //opserr << " FE f_pred = " << f_pred << endln; 03192 03193 03194 03195 stresstensor intersection_stress = start_stress; // added 20april2000 for forward euler 03196 03197 stresstensor elpl_start_stress = start_stress; 03198 03199 stresstensor true_stress_increment = stress_increment; 03200 03201 straintensor El_strain_increment; 03202 03203 03204 03205 if ( f_start <= 0 && f_pred <= 0 || f_start > f_pred ) 03206 03207 { 03208 03209 //Updating elastic strain increment 03210 03211 straintensor estrain = ElasticPredictorEPS.getElasticStrain(); 03212 03213 straintensor tstrain = ElasticPredictorEPS.getStrain(); 03214 03215 estrain = estrain + strain_incr; 03216 03217 tstrain = tstrain + strain_incr; 03218 03219 ElasticPredictorEPS.setElasticStrain( estrain ); 03220 03221 ElasticPredictorEPS.setStrain( tstrain ); 03222 03223 ElasticPredictorEPS.setdElasticStrain( strain_incr ); 03224 03225 03226 03227 //Evolve parameters like void ratio (e) according to elastic strain and elastic stress--for MD model specifically 03228 03229 //double Delta_lambda = 0.0; 03230 03231 //material_point.EL->UpdateVar( &ElasticPredictorEPS, 1); 03232 03233 // Update E_Young and e according to current stress state before evaluate ElasticStiffnessTensor 03234 03235 if ( getELT1() ) { 03236 03237 //getELT1()->updateEeDm(&ElasticPredictorEPS, st_vol, 0.0); 03238 03239 getELT1()->updateEeDm(&ElasticPredictorEPS, -st_vol, 0.0); 03240 03241 } 03242 03243 03244 03245 //opserr <<" strain_increment.Iinvariant1() " << strain_increment.Iinvariant1() << endlnn; 03246 03247 ElasticPredictorEPS.setEep(E); 03248 03249 return ElasticPredictorEPS; 03250 03251 } 03252 03253 03254 03255 if ( f_start <= 0 && f_pred > 0 ) 03256 03257 { 03258 03259 intersection_stress = 03260 03261 yield_surface_cross( start_stress, elastic_predictor_stress); 03262 03263 //opserr << " start_stress: " << start_stress << endlnn; 03264 03265 //opserr << " Intersection_stress: " << intersection_stress << endlnn; 03266 03267 03268 03269 IntersectionEPS.setStress( intersection_stress ); 03270 03271 //intersection_stress.reportshort("Intersection stress \n"); 03272 03273 03274 03275 elpl_start_stress = intersection_stress; 03276 03277 //elpl_start_stress.reportshortpqtheta("elpl start stress AFTER \n"); 03278 03279 03280 03281 true_stress_increment = elastic_predictor_stress - elpl_start_stress; 03282 03283 //true_stress_increment.null_indices(); 03284 03285 03286 03287 stresstensor EstressIncr = intersection_stress- start_stress; 03288 03289 03290 03291 //forwardEPS.setStress( elpl_start_stress ); 03292 03293 //Should only count on that elastic portion, not st_vol... 03294 03295 if ( getELT1() ) { 03296 03297 El_strain_increment = D("ijpq") * EstressIncr("pq"); 03298 03299 double st_vol_El_incr = El_strain_increment.Iinvariant1(); 03300 03301 03302 03303 //opserr << " FE crossing update... "; 03304 03305 getELT1()->updateEeDm(&IntersectionEPS, -st_vol_El_incr, 0.0); 03306 03307 //opserr << " FE crossing update... "; 03308 03309 getELT1()->updateEeDm(&forwardEPS, -st_vol_El_incr, 0.0); 03310 03311 } 03312 03313 03314 03315 } 03316 03317 03318 03319 03320 03321 //forwardEPS.setStress( elpl_start_stress ); 03322 03323 03324 03325 //opserr <<"elpl start eps: " << forwardEPS; 03326 03327 //double f_cross = this->getYS()->f( &forwardEPS ); 03328 03329 //opserr << " ####### f_cross = " << f_cross << "\n"; 03330 03331 03332 03333 //set the initial value of D once the current stress hits the y.s. for Manzari-Dafalias Model 03334 03335 //if ( f_start <= 0 && f_pred > 0 ) 03336 03337 // material_point.EL->setInitD(&forwardEPS); 03338 03339 //opserr << " inside ConstitutiveDriver after setInitD " << forwardEPS; 03340 03341 03342 03343 03344 03345 // pulling out some tensor and double definitions 03346 03347 //tensor dFods( 2, def_dim_2, 0.0); 03348 03349 //tensor dQods( 2, def_dim_2, 0.0); 03350 03351 stresstensor dFods; 03352 03353 stresstensor dQods; 03354 03355 // stresstensor s; // deviator 03356 03357 tensor H( 2, def_dim_2, 0.0); 03358 03359 tensor temp1( 2, def_dim_2, 0.0); 03360 03361 tensor temp2( 2, def_dim_2, 0.0); 03362 03363 double lower = 0.0; 03364 03365 tensor temp3( 2, def_dim_2, 0.0); 03366 03367 03368 03369 double Delta_lambda = 0.0; 03370 03371 double h_s[4] = {0.0, 0.0, 0.0, 0.0}; 03372 03373 double xi_s[4] = {0.0, 0.0, 0.0, 0.0}; 03374 03375 stresstensor h_t[4]; 03376 03377 stresstensor xi_t[4]; 03378 03379 double hardMod_ = 0.0; 03380 03381 03382 03383 //double Dq_ast = 0.0; 03384 03385 //double q_ast_entry = 0.0; 03386 03387 //double q_ast = 0.0; 03388 03389 03390 03391 stresstensor plastic_stress; 03392 03393 straintensor plastic_strain; 03394 03395 stresstensor elastic_plastic_stress; 03396 03397 // ::printf("\n...... felpred = %lf\n",felpred); 03398 03399 03400 03401 if ( f_pred >= 0 ) { 03402 03403 03404 03405 03406 03407 //dFods = getYS()->dFods( &forwardEPS ); 03408 03409 //dQods = getPS()->dQods( &forwardEPS ); 03410 03411 dFods = getYS()->dFods( &IntersectionEPS ); 03412 03413 dQods = getPS()->dQods( &IntersectionEPS ); 03414 03415 03416 03417 //opserr << "dF/ds" << dFods << endlnn; 03418 03419 //opserr << "dQ/ds" << dQods << endlnn; 03420 03421 03422 03423 // Tensor H_kl ( eq. 5.209 ) W.F. Chen 03424 03425 H = E("ijkl")*dQods("kl"); //E_ijkl * R_kl 03426 03427 H.null_indices(); 03428 03429 temp1 = dFods("ij") * E("ijkl"); // L_ij * E_ijkl 03430 03431 temp1.null_indices(); 03432 03433 temp2 = temp1("ij")*dQods("ij"); // L_ij * E_ijkl * R_kl 03434 03435 temp2.null_indices(); 03436 03437 lower = temp2.trace(); 03438 03439 03440 03441 // Evaluating the hardening modulus: sum of (df/dq*) * qbar 03442 03443 03444 03445 hardMod_ = 0.0; 03446 03447 //Of 1st scalar internal vars 03448 03449 if ( getELS1() ) { 03450 03451 h_s[0] = getELS1()->h_s(&IntersectionEPS, getPS()); 03452 03453 xi_s[0] = getYS()->xi_s1( &IntersectionEPS ); 03454 03455 hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 03456 03457 } 03458 03459 03460 03461 //Of 2nd scalar internal vars 03462 03463 if ( getELS2() ) { 03464 03465 h_s[1] = getELS2()->h_s( &IntersectionEPS, getPS()); 03466 03467 xi_s[1] = getYS()->xi_s2( &IntersectionEPS ); 03468 03469 hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 03470 03471 } 03472 03473 03474 03475 //Of 3rd scalar internal vars 03476 03477 if ( getELS3() ) { 03478 03479 h_s[2] = getELS3()->h_s( &IntersectionEPS, getPS()); 03480 03481 xi_s[2] = getYS()->xi_s3( &IntersectionEPS ); 03482 03483 hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 03484 03485 } 03486 03487 03488 03489 //Of 4th scalar internal vars 03490 03491 if ( getELS4() ) { 03492 03493 h_s[3] = getELS4()->h_s( &IntersectionEPS, getPS()); 03494 03495 xi_s[3] = getYS()->xi_s4( &IntersectionEPS ); 03496 03497 hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 03498 03499 } 03500 03501 03502 03503 //Of tensorial internal var 03504 03505 // 1st tensorial var 03506 03507 if ( getELT1() ) { 03508 03509 h_t[0] = getELT1()->h_t(&IntersectionEPS, getPS()); 03510 03511 xi_t[0] = getYS()->xi_t1( &IntersectionEPS ); 03512 03513 tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 03514 03515 hardMod_ = hardMod_ + hm.trace(); 03516 03517 } 03518 03519 03520 03521 // 2nd tensorial var 03522 03523 if ( getELT2() ) { 03524 03525 h_t[1] = getELT2()->h_t( &IntersectionEPS, getPS()); 03526 03527 xi_t[1] = getYS()->xi_t2( &IntersectionEPS ); 03528 03529 tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 03530 03531 hardMod_ = hardMod_ + hm.trace(); 03532 03533 } 03534 03535 03536 03537 // 3rd tensorial var 03538 03539 if ( getELT3() ) { 03540 03541 h_t[2] = getELT3()->h_t( &IntersectionEPS, getPS()); 03542 03543 xi_t[2] = getYS()->xi_t3( &IntersectionEPS ); 03544 03545 tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 03546 03547 hardMod_ = hardMod_ + hm.trace(); 03548 03549 } 03550 03551 03552 03553 // 4th tensorial var 03554 03555 if ( getELT4() ) { 03556 03557 h_t[3] = getELT4()->h_t(&IntersectionEPS, getPS()); 03558 03559 xi_t[3] = getYS()->xi_t4( &IntersectionEPS ); 03560 03561 tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 03562 03563 hardMod_ = hardMod_ + hm.trace(); 03564 03565 } 03566 03567 03568 03569 // Subtract accumulated hardMod_ from lower 03570 03571 lower = lower - hardMod_; 03572 03573 //opserr << "FE : lower " << lower << " hardMod_ " << hardMod_ << endln; 03574 03575 //Calculating Kp according to Kp = - (df/dq*) * qbar 03576 03577 //double Kp = material_point.EL->getKp(&forwardEPS, norm_dQods); 03578 03579 //Kp = 0.0; 03580 03581 //opserr << endlnn << ">>>>>>>>> Lower = " << lower << endlnn; 03582 03583 //lower = lower + Kp; 03584 03585 //opserr << endlnn << ">>>>>>>>> Kp = " << Kp << endlnn; 03586 03587 03588 03589 //opserr << " stress_increment "<< stress_increment << endlnn; 03590 03591 //opserr << " true_stress_increment "<< true_stress_increment << endlnn; 03592 03593 03594 03595 temp3 = dFods("ij") * true_stress_increment("ij"); // L_ij * E_ijkl * d e_kl (true ep strain increment) 03596 03597 temp3.null_indices(); 03598 03599 //opserr << " temp3.trace() -- true_stress_incr " << temp3.trace() << endln; 03600 03601 //GZ temp3 = temp1("ij")*strain_incr("ij"); 03602 03603 //GZ temp3.null_indices(); 03604 03605 //opserr << " temp3.trace() " << temp3.trace() << endlnn; 03606 03607 Delta_lambda = (temp3.trace())/lower; 03608 03609 //opserr << "FE: Delta_lambda " << Delta_lambda << endln; 03610 03611 if (Delta_lambda<0.0) Delta_lambda=0.0; 03612 03613 03614 03615 plastic_stress = H("kl") * Delta_lambda; 03616 03617 plastic_strain = dQods("kl") * Delta_lambda; // plastic strain increment 03618 03619 plastic_stress.null_indices(); 03620 03621 plastic_strain.null_indices(); 03622 03623 //opserr << " Delta_lambda " << Delta_lambda << "plastic_stress = " << plastic_stress << endln; 03624 03625 //opserr << "plastic_stress = " << plastic_stress << endlnn; 03626 03627 //opserr << "plastic_strain = " << plastic_strain << endlnn; 03628 03629 //opserr << "plastic_strain I1= " << plastic_strain.Iinvariant1() << endlnn; 03630 03631 //opserr << "plastic_strain vol " << Delta_lambda * ( forwardEPS.getScalarVar( 2 ) )<< endlnn ; 03632 03633 //opserr << " q=" << Delta_lambda * dQods.q_deviatoric()<< endlnn; 03634 03635 //plastic_stress.reportshort("plastic stress (with delta_lambda)\n"); 03636 03637 03638 03639 elastic_plastic_stress = elastic_predictor_stress - plastic_stress; 03640 03641 //elastic_plastic_stress.reportshortpqtheta("FE elastic plastic stress \n"); 03642 03643 03644 03645 //calculating elatic strain increment 03646 03647 //stresstensor dstress_el = elastic_plastic_stress - start_stress; 03648 03649 //straintensor elastic_strain = D("ijpq") * dstress_el("pq"); 03650 03651 straintensor elastic_strain = strain_incr - plastic_strain; // elastic strain increment 03652 03653 //opserr << "elastic_strain I1=" << elastic_strain.Iinvariant1() << endlnn; 03654 03655 //opserr << "elastic_strain " << elastic_strain << endlnn; 03656 03657 //opserr << "strain increment I1=" << strain_increment.Iinvariant1() << endlnn; 03658 03659 //opserr << "strain increment " << strain_increment << endlnn; 03660 03661 03662 03663 straintensor estrain = forwardEPS.getElasticStrain(); //get old elastic strain 03664 03665 straintensor pstrain = forwardEPS.getPlasticStrain(); //get old plastic strain 03666 03667 //straintensor pstrain = forwardEPS.getStrain_commit() - estrain; //get old plastic strain 03668 03669 03670 03671 straintensor tstrain = forwardEPS.getStrain(); //get old total strain 03672 03673 pstrain = pstrain + plastic_strain; 03674 03675 estrain = estrain + elastic_strain; 03676 03677 tstrain = tstrain + elastic_strain + plastic_strain; 03678 03679 03680 03681 //Setting de_p, de_e, total plastic, elastic strain, and total strain 03682 03683 forwardEPS.setdPlasticStrain( plastic_strain ); 03684 03685 forwardEPS.setdElasticStrain( elastic_strain ); 03686 03687 forwardEPS.setPlasticStrain( pstrain ); 03688 03689 forwardEPS.setElasticStrain( estrain ); 03690 03691 forwardEPS.setStrain( tstrain ); 03692 03693 03694 03695 //================================================================ 03696 03697 //Generating Eep using dQods at the intersection point 03698 03699 dFods = getYS()->dFods( &IntersectionEPS ); 03700 03701 dQods = getPS()->dQods( &IntersectionEPS ); 03702 03703 03704 03705 tensor upperE1 = E("pqkl")*dQods("kl"); 03706 03707 upperE1.null_indices(); 03708 03709 tensor upperE2 = dFods("ij")*E("ijmn"); 03710 03711 upperE2.null_indices(); 03712 03713 03714 03715 //tensor upperE = upperE1("pq") * upperE1("mn"); // Bug found, Zhao Cheng Jan13, 2004 03716 03717 tensor upperE = upperE1("pq") * upperE2("mn"); 03718 03719 upperE.null_indices(); 03720 03721 03722 03723 /*//temp2 = upperE2("ij")*dQods("ij"); // L_ij * E_ijkl * R_kl 03724 03725 temp2.null_indices(); 03726 03727 lower = temp2.trace(); 03728 03729 03730 03731 03732 03733 // Evaluating the hardening modulus: sum of (df/dq*) * qbar 03734 03735 03736 03737 hardMod_ = 0.0; 03738 03739 //Of 1st scalar internal vars 03740 03741 if ( getELS1() ) { 03742 03743 h_s[0] = getELS1()->h_s( &IntersectionEPS, getPS()); 03744 03745 xi_s[0] = getYS()->xi_s1( &IntersectionEPS ); 03746 03747 hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 03748 03749 } 03750 03751 03752 03753 //Of 2nd scalar internal vars 03754 03755 if ( getELS2() ) { 03756 03757 h_s[1] = getELS2()->h_s( &IntersectionEPS, getPS()); 03758 03759 xi_s[1] = getYS()->xi_s2( &IntersectionEPS ); 03760 03761 hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 03762 03763 } 03764 03765 03766 03767 //Of 3rd scalar internal vars 03768 03769 if ( getELS3() ) { 03770 03771 h_s[2] = getELS3()->h_s( &IntersectionEPS, getPS()); 03772 03773 xi_s[2] = getYS()->xi_s3( &IntersectionEPS ); 03774 03775 hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 03776 03777 } 03778 03779 03780 03781 //Of 4th scalar internal vars 03782 03783 if ( getELS4() ) { 03784 03785 h_s[3] = getELS4()->h_s( &IntersectionEPS, getPS()); 03786 03787 xi_s[3] = getYS()->xi_s4( &IntersectionEPS ); 03788 03789 hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 03790 03791 } 03792 03793 03794 03795 //Of tensorial internal var 03796 03797 // 1st tensorial var 03798 03799 if ( getELT1() ) { 03800 03801 h_t[0] = getELT1()->h_t(&IntersectionEPS, getPS()); 03802 03803 xi_t[0] = getYS()->xi_t1( &IntersectionEPS ); 03804 03805 tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 03806 03807 hardMod_ = hardMod_ + hm.trace(); 03808 03809 } 03810 03811 03812 03813 // 2nd tensorial var 03814 03815 if ( getELT2() ) { 03816 03817 h_t[1] = getELT2()->h_t( &IntersectionEPS, getPS()); 03818 03819 xi_t[1] = getYS()->xi_t2( &IntersectionEPS ); 03820 03821 tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 03822 03823 hardMod_ = hardMod_ + hm.trace(); 03824 03825 } 03826 03827 03828 03829 // 3rd tensorial var 03830 03831 if ( getELT3() ) { 03832 03833 h_t[2] = getELT3()->h_t(&IntersectionEPS, getPS()); 03834 03835 xi_t[2] = getYS()->xi_t3( &IntersectionEPS ); 03836 03837 tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 03838 03839 hardMod_ = hardMod_ + hm.trace(); 03840 03841 } 03842 03843 03844 03845 // 4th tensorial var 03846 03847 if ( getELT4() ) { 03848 03849 h_t[3] = getELT4()->h_t( &IntersectionEPS, getPS()); 03850 03851 xi_t[3] = getYS()->xi_t4( &IntersectionEPS ); 03852 03853 tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 03854 03855 hardMod_ = hardMod_ + hm.trace(); 03856 03857 } 03858 03859 03860 03861 // Subtract accumulated hardMod_ from lower 03862 03863 03864 03865 lower = lower - hardMod_; 03866 03867 */ 03868 03869 03870 03871 tensor Ep = upperE*(1./lower); 03872 03873 03874 03875 // elastoplastic constitutive tensor 03876 03877 double h_L = 0.0; // Bug fixed Joey 07-21-02 added h(L) function 03878 03879 if ( Delta_lambda > 0 ) h_L = 1.0; 03880 03881 //opserr << " h_L = " << h_L << "\n"; 03882 03883 //Eep = Eep - Ep*h_L; // Bug found, Zhao Cheng Jan13, 2004 03884 03885 Eep = E - Ep*h_L; 03886 03887 03888 03889 //opserr <<" after calculation---Eep.rank()= " << Eep.rank() <<endlnn; 03890 03891 //Eep.printshort(" IN template "); 03892 03893 03894 03895 //--// before the surface is been updated ! 03896 03897 //--// f_Final = Criterion.f(elastic_plastic_stress); 03898 03899 //-- 03900 03901 //-- q_ast_entry = Criterion.kappa_get(elastic_plastic_stress); 03902 03903 //-- 03904 03905 //--// h_ = h(elastic_plastic_stress); 03906 03907 //-- Dq_ast = Delta_lambda * h_ * just_this_PP; 03908 03909 //--// if (Dq_ast < 0.0 ) Dq_ast = 0.0; 03910 03911 //--// Dq_ast = fabs(Delta_lambda * h_ * just_this_PP); // because of softening response... 03912 03913 //--//::printf(" h_=%.6e q_ast_entry=%.6e Dq_ast=%.6e Delta_lambda=%.6e\n", h_, q_ast_entry, Dq_ast, Delta_lambda); 03914 03915 //-- 03916 03917 //-- current_lambda_set(Delta_lambda); 03918 03919 //-- 03920 03921 //-- q_ast = q_ast_entry + Dq_ast; 03922 03923 //-- Criterion.kappa_set( elastic_plastic_stress, q_ast); 03924 03925 //--//::fprintf(stdout," Criterion.kappa_get(elastic_plastic_stress)=%.8e\n",Criterion.kappa_get(elastic_plastic_stress)); 03926 03927 //--//::fprintf(stderr," Criterion.kappa_get(elastic_plastic_stress)=%.8e\n",Criterion.kappa_get(elastic_plastic_stress)); 03928 03929 //-- 03930 03931 //-- 03932 03933 //--//::fprintf(stdout," ######## predictor --> q_ast_entry=%.8e Dq_ast=%.8e q_ast=%.8e\n",q_ast_entry, Dq_ast, q_ast); 03934 03935 //--//::fprintf(stderr," ######## predictor --> q_ast_entry=%.8e Dq_ast=%.8e q_ast=%.8e\n",q_ast_entry, Dq_ast, q_ast); 03936 03937 //-- 03938 03939 //--//::fprintf(stdout,"ForwardEulerStress IN Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress)); 03940 03941 //--//::fprintf(stderr,"ForwardEulerStress IN Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress)); 03942 03943 //-- 03944 03945 03946 03947 //new EPState in terms of stress 03948 03949 forwardEPS.setStress(elastic_plastic_stress); 03950 03951 //opserr <<"inside constitutive driver: forwardEPS "<< forwardEPS; 03952 03953 03954 03955 forwardEPS.setEep(Eep); 03956 03957 //forwardEPS.getEep().printshort(" after set"); 03958 03959 03960 03961 //Before update all the internal vars 03962 03963 double f_forward = getYS()->f( &forwardEPS ); 03964 03965 //opserr << "\n************ Before f_forward = " << f_forward << "\n"; 03966 03967 03968 03969 //Evolve the surfaces and hardening vars 03970 03971 int NS = forwardEPS.getNScalarVar(); 03972 03973 int NT = forwardEPS.getNTensorVar(); 03974 03975 03976 03977 double dS= 0.0; 03978 03979 double S = 0.0; 03980 03981 03982 03983 int ii; 03984 03985 for (ii = 1; ii <= NS; ii++) { 03986 03987 dS = Delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var 03988 03989 S = forwardEPS.getScalarVar(ii); //Guanzhou Mar2005 // Get the old value of the scalar internal var 03990 03991 forwardEPS.setScalarVar(ii, S + dS ); // Update internal scalar var 03992 03993 } 03994 03995 03996 03997 stresstensor dT; 03998 03999 stresstensor T; 04000 04001 stresstensor new_T; 04002 04003 04004 04005 for (ii = 1; ii <= NT; ii++) { 04006 04007 dT = h_t[ii-1]*Delta_lambda ; // Increment to the tensor internal var 04008 04009 T = forwardEPS.getTensorVar(ii); //Guanzhou Mar2005 // Get the old value of the tensor internal var 04010 04011 new_T = T + dT; 04012 04013 forwardEPS.setTensorVar(ii, new_T ); 04014 04015 } 04016 04017 04018 04019 // Update E_Young and e according to current stress state 04020 04021 //straintensor pl_strain_increment; 04022 04023 int err = 0; 04024 04025 if ( getELT1() ) { 04026 04027 //pl_strain_increment = strain_increment - El_strain_increment; 04028 04029 double st_vol_pl = plastic_strain.Iinvariant1(); 04030 04031 //err = getELT1()->updateEeDm(&forwardEPS, st_vol_pl, Delta_lambda); 04032 04033 //opserr << "pl_vol= " << st_vol_pl << "|update before FE \n"; 04034 04035 //opserr << " FE Pl update..."; 04036 04037 // D > 0 compressive -> Iinv > 0 -> de < 0 correct! 04038 04039 err = getELT1()->updateEeDm(&forwardEPS, st_vol_pl, Delta_lambda); 04040 04041 } 04042 04043 04044 04045 //tensor tempx = plastic_strain("ij") * plastic_strain("ij"); 04046 04047 //double tempxd = tempx.trace(); 04048 04049 //double e_eq = pow( 2.0 * tempxd / 3.0, 0.5 ); 04050 04052 04053 // 04054 04055 //double dalfa1 = e_eq * 10; 04056 04057 //double alfa1 = forwardEPS.getScalarVar(1); 04058 04059 04060 04061 04062 04063 04064 04065 //opserr << "UpdateAllVars " << forwardEPS<< endlnn; 04066 04067 04068 04069 //After update all the internal vars 04070 04071 f_forward = getYS()->f( &forwardEPS ); 04072 04073 //opserr << "\n************ After f_forward = " << f_forward << "\n"; 04074 04075 04076 04077 04078 04079 } 04080 04081 04082 04083 //::fprintf(stderr,"ForwardEulerStress EXIT Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress)); 04084 04085 //forwardEPS.setConverged(TRUE); 04086 04087 04088 04089 //double p = (forwardEPS.getStress()).p_hydrostatic(); 04090 04091 //double ec = (forwardEPS.getec()) - (forwardEPS.getLam()) * log( p / (forwardEPS.getpo()) ); 04092 04093 //double st = (forwardEPS.getStrain()).Iinvariant1(); 04094 04095 //double pl_s = (forwardEPS.getPlasticStrain()).Iinvariant1(); 04096 04097 //double dpl_s = (forwardEPS.getdPlasticStrain()).Iinvariant1(); 04098 04099 //opserr << "P FE p=" << p << " ec " << ec << " e " << forwardEPS.gete() << " psi " << (forwardEPS.gete() - ec) << " strain " << st << " t_pl " << pl_s << " d_pl " << dpl_s << "\n"; 04100 04101 forwardEPS.Delta_lambda = Delta_lambda; 04102 04103 04104 04105 return forwardEPS; 04106 04107 } 04108 04109 04110 04111 //================================================================================ 04112 04113 // Starting EPState using Semi Backward Euler Starting Point 04114 04115 //================================================================================ 04116 04117 EPState Template3Dep::SemiBackwardEulerEPState( const straintensor &strain_increment) 04118 04119 { 04120 04121 stresstensor start_stress; 04122 04123 EPState SemibackwardEPS( *(this->getEPS()) ); 04124 04125 start_stress = SemibackwardEPS.getStress(); 04126 04127 04128 04129 // building elasticity tensor 04130 04131 //tensor E = Criterion.StiffnessTensorE(Ey,nu); 04132 04133 tensor E = ElasticStiffnessTensor(); 04134 04135 // building compliance tensor 04136 04137 // tensor D = Criterion.ComplianceTensorD(Ey,nu); 04138 04139 04140 04141 // pulling out some tensor and double definitions 04142 04143 tensor dFods(2, def_dim_2, 0.0); 04144 04145 tensor dQods(2, def_dim_2, 0.0); 04146 04147 tensor temp2(2, def_dim_2, 0.0); 04148 04149 double lower = 0.0; 04150 04151 double Delta_lambda = 0.0; 04152 04153 04154 04155 EPState E_Pred_EPS( *(this->getEPS()) ); 04156 04157 04158 04159 straintensor strain_incr = strain_increment; 04160 04161 stresstensor stress_increment = E("ijkl")*strain_incr("kl"); 04162 04163 stress_increment.null_indices(); 04164 04165 // stress_increment.reportshort("stress Increment\n"); 04166 04167 04168 04169 04170 04171 stresstensor plastic_stress; 04172 04173 stresstensor elastic_predictor_stress; 04174 04175 stresstensor elastic_plastic_stress; 04176 04177 //.. double Felplpredictor = 0.0; 04178 04179 04180 04181 double h_s[4] = {0.0, 0.0, 0.0, 0.0}; 04182 04183 double xi_s[4] = {0.0, 0.0, 0.0, 0.0}; 04184 04185 stresstensor h_t[4]; 04186 04187 stresstensor xi_t[4]; 04188 04189 double hardMod_ = 0.0; 04190 04191 04192 04193 double S = 0.0; 04194 04195 double dS = 0.0; 04196 04197 stresstensor T; 04198 04199 stresstensor dT; 04200 04201 //double Dq_ast = 0.0; 04202 04203 //double q_ast_entry = 0.0; 04204 04205 //double q_ast = 0.0; 04206 04207 04208 04209 elastic_predictor_stress = start_stress + stress_increment; 04210 04211 //.. elastic_predictor_stress.reportshort("ELASTIC PREDICTOR stress\n"); 04212 04213 E_Pred_EPS.setStress( elastic_predictor_stress ); 04214 04215 04216 04217 // double f_start = Criterion.f(start_stress); 04218 04219 // ::printf("SEMI BE############## f_start = %.10e\n",f_start); 04220 04221 double f_pred = getYS()->f( &E_Pred_EPS ); 04222 04223 //::printf("SEMI BE############## f_pred = %.10e\n",f_pred); 04224 04225 04226 04227 // second of alternative predictors as seen in MAC page 170-171 04228 04229 if ( f_pred >= 0.0 ) 04230 04231 { 04232 04233 //el_or_pl_range(1); // set to 1 ( plastic range ) 04234 04235 // PREDICTOR FASE 04236 04237 //.. ::printf("\n\npredictor part step_counter = %d\n\n", step_counter); 04238 04239 04240 04241 dFods = getYS()->dFods( &E_Pred_EPS ); 04242 04243 dQods = getPS()->dQods( &E_Pred_EPS ); 04244 04245 //.... dFods.print("a","dF/ds "); 04246 04247 //.... dQods.print("a","dQ/ds "); 04248 04249 04250 04251 temp2 = (dFods("ij")*E("ijkl"))*dQods("kl"); 04252 04253 temp2.null_indices(); 04254 04255 lower = temp2.trace(); 04256 04257 04258 04259 // Evaluating the hardening modulus: sum of (df/dq*) * qbar 04260 04261 04262 04263 //Of scalar internal var 04264 04265 hardMod_ = 0.0; 04266 04267 04268 04269 //Of 1st scalar internal vars 04270 04271 if ( getELS1() ) { 04272 04273 h_s[0] = getELS1()->h_s( &E_Pred_EPS, getPS()); 04274 04275 xi_s[0] = getYS()->xi_s1( &E_Pred_EPS ); 04276 04277 hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 04278 04279 } 04280 04281 04282 04283 //Of 2nd scalar internal vars 04284 04285 if ( getELS2() ) { 04286 04287 h_s[1] = getELS2()->h_s( &E_Pred_EPS, getPS()); 04288 04289 xi_s[1] = getYS()->xi_s2( &E_Pred_EPS ); 04290 04291 hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 04292 04293 } 04294 04295 04296 04297 //Of 3rd scalar internal vars 04298 04299 if ( getELS3() ) { 04300 04301 h_s[2] = getELS3()->h_s(&E_Pred_EPS, getPS()); 04302 04303 xi_s[2] = getYS()->xi_s3( &E_Pred_EPS ); 04304 04305 hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 04306 04307 } 04308 04309 04310 04311 //Of 4th scalar internal vars 04312 04313 if ( getELS4() ) { 04314 04315 h_s[3] = getELS4()->h_s( &E_Pred_EPS, getPS()); 04316 04317 xi_s[3] = getYS()->xi_s4( &E_Pred_EPS ); 04318 04319 hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 04320 04321 } 04322 04323 04324 04325 //Of tensorial internal var 04326 04327 // 1st tensorial var 04328 04329 if ( getELT1() ) { 04330 04331 h_t[0] = getELT1()->h_t(&E_Pred_EPS, getPS()); 04332 04333 xi_t[0] = getYS()->xi_t1( &E_Pred_EPS ); 04334 04335 tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 04336 04337 hardMod_ = hardMod_ + hm.trace(); 04338 04339 } 04340 04341 04342 04343 // 2nd tensorial var 04344 04345 if ( getELT2() ) { 04346 04347 h_t[1] = getELT2()->h_t(&E_Pred_EPS, getPS()); 04348 04349 xi_t[1] = getYS()->xi_t2( &E_Pred_EPS ); 04350 04351 tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 04352 04353 hardMod_ = hardMod_ + hm.trace(); 04354 04355 } 04356 04357 04358 04359 // 3rd tensorial var 04360 04361 if ( getELT3() ) { 04362 04363 h_t[2] = getELT3()->h_t(&E_Pred_EPS, getPS()); 04364 04365 xi_t[2] = getYS()->xi_t3( &E_Pred_EPS ); 04366 04367 tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 04368 04369 hardMod_ = hardMod_ + hm.trace(); 04370 04371 } 04372 04373 04374 04375 // 4th tensorial var 04376 04377 if ( getELT4() ) { 04378 04379 h_t[3] = getELT4()->h_t(&E_Pred_EPS, getPS()); 04380 04381 xi_t[3] = getYS()->xi_t4( &E_Pred_EPS ); 04382 04383 tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 04384 04385 hardMod_ = hardMod_ + hm.trace(); 04386 04387 } 04388 04389 04390 04391 // Subtract accumulated hardMod_ from lower 04392 04393 lower = lower - hardMod_; 04394 04395 04396 04397 //h_s = material_point.ELS1->h_s( &E_Pred_EPS, material_point.PS ); 04398 04399 //xi_s = material_point.YS->xi_s1( &E_Pred_EPS ); 04400 04401 //hardMod_ = h_s * xi_s; 04402 04403 //lower = lower - hardMod_; 04404 04405 04406 04408 04409 //h_t = material_point.ELT1->h_t(&E_Pred_EPS, material_point.PS); 04410 04411 //xi_t = material_point.YS->xi_t1( &E_Pred_EPS ); 04412 04413 //tensor hm = h_t("ij") * xi_t("ij"); 04414 04415 //hardMod_ = hm.trace(); 04416 04417 //lower = lower - hardMod_; 04418 04419 04420 04421 04422 04423 Delta_lambda = f_pred/lower; 04424 04425 if ( Delta_lambda < 0.0 ) 04426 04427 { 04428 04429 //::fprintf(stderr,"\nP!\n"); 04430 04431 } 04432 04433 plastic_stress = (E("ijkl")*dQods("kl"))*(-Delta_lambda); 04434 04435 plastic_stress.null_indices(); 04436 04437 //.. plastic_stress.reportshort("plastic stress predictor II\n"); 04438 04439 //.. elastic_predictor_stress.reportshort("elastic predictor stress \n"); 04440 04441 elastic_plastic_stress = elastic_predictor_stress + plastic_stress; 04442 04443 elastic_plastic_stress.null_indices(); 04444 04445 04446 04447 SemibackwardEPS.setStress( elastic_plastic_stress ); 04448 04449 04450 04452 04453 //S = SemibackwardEPS.getScalarVar(1); // Get the old value of the internal var 04454 04455 //h_s = material_point.ELS1->h_s( &SemibackwardEPS, material_point.PS ); 04456 04457 //dS = Delta_lambda * h_s ; // Increment to the internal scalar var 04458 04459 //h_t = material_point.ELT1->h_t( &SemibackwardEPS, material_point.PS ); 04460 04461 //dT = Delta_lambda * h_t ; // Increment to the internal tensorial var 04462 04463 04464 04465 //Evolve the surfaces and hardening vars 04466 04467 int NS = SemibackwardEPS.getNScalarVar(); 04468 04469 int NT = SemibackwardEPS.getNTensorVar(); 04470 04471 04472 04473 int ii; 04474 04475 for (ii = 1; ii <= NS; ii++) { 04476 04477 dS = Delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var 04478 04479 S = SemibackwardEPS.getScalarVar(ii); // Get the old value of the scalar internal var 04480 04481 SemibackwardEPS.setScalarVar(ii, S + dS ); // Update internal scalar var 04482 04483 } 04484 04485 04486 04487 04488 04490 04492 04494 04495 //SemibackwardEPS.setScalarVar(1, S + dS ); 04496 04497 04498 04499 //stresstensor new_T = T + dT; 04500 04501 //SemibackwardEPS.setTensorVar(1, new_T ); 04502 04503 04504 04505 stresstensor new_T; 04506 04507 04508 04509 for (ii = 1; ii <= NT; ii++) { 04510 04511 dT = h_t[ii-1] * Delta_lambda; // Increment to the tensor internal var 04512 04513 T = SemibackwardEPS.getTensorVar(ii); // Get the old value of the tensor internal var 04514 04515 new_T = T + dT; 04516 04517 SemibackwardEPS.setTensorVar(ii, new_T ); 04518 04519 } 04520 04521 04522 04523 04524 04525 //return elastic_plastic_stress; 04526 04527 return SemibackwardEPS; 04528 04529 } 04530 04531 return E_Pred_EPS; 04532 04533 } 04534 04535 04536 04537 04538 04539 04540 04541 04542 04543 04544 04545 //================================================================================ 04546 04547 // New EPState using Backward Euler Algorithm 04548 04549 //================================================================================ 04550 04551 EPState Template3Dep::BackwardEulerEPState( const straintensor &strain_increment) 04552 04553 04554 04555 { 04556 04558 04559 04560 04561 // Volumetric strain 04562 04563 //double st_vol = strain_increment.p_hydrostatic(); 04564 04565 double st_vol = strain_increment.Iinvariant1(); //- = compressive 04566 04567 04568 04569 //EPState to be returned, it can be elastic or elastic-plastic EPState 04570 04571 EPState backwardEPS( * (this->getEPS()) ); 04572 04573 EPState startEPS( *(this->getEPS()) ); 04574 04575 04576 04577 stresstensor start_stress = startEPS.getStress(); //Guanzhou Mar2005 04578 04579 04580 04581 //Output for plotting 04582 04583 //opserr.precision(5); 04584 04585 //opserr.width(10); 04586 04587 //opserr << " strain_increment " << strain_increment << "\n"; 04588 04589 04590 04591 //opserr.precision(5); 04592 04593 //opserr.width(10); 04594 04595 //opserr << "start_stress " << start_stress; 04596 04597 04598 04599 // Pulling out some tensor and double definitions 04600 04601 tensor I2("I", 2, def_dim_2); 04602 04603 tensor I_ijkl("I", 4, def_dim_4); 04604 04605 I_ijkl = I2("ij")*I2("kl"); 04606 04607 I_ijkl.null_indices(); 04608 04609 tensor I_ikjl("I", 4, def_dim_4); 04610 04611 I_ikjl = I_ijkl.transpose0110(); 04612 04613 04614 04615 //double Ey = Criterion.E(); 04616 04617 //double nu = Criterion.nu(); 04618 04619 //tensor E = StiffnessTensorE(Ey,nu); 04620 04621 tensor E = ElasticStiffnessTensor(); 04622 04623 // tensor D = ComplianceTensorD(Ey,nu); 04624 04625 // stresstensor REAL_start_stress = start_stress; 04626 04627 04628 04629 tensor dFods(2, def_dim_2, 0.0); 04630 04631 tensor dQods(2, def_dim_2, 0.0); 04632 04633 // tensor dQodsextended(2, def_dim_2, 0.0); 04634 04635 // tensor d2Qodqast(2, def_dim_2, 0.0); 04636 04637 tensor temp2(2, def_dim_2, 0.0); 04638 04639 double lower = 0.0; 04640 04641 double Delta_lambda = 0.0; // Lambda 04642 04643 double delta_lambda = 0.0; // dLambda 04644 04645 04646 04647 double Felplpredictor = 0.0; 04648 04649 //Kai double absFelplpredictor = 0.0; 04650 04651 // double Ftolerance = pow(d_macheps(),(1.0/2.0))*1000000.00; //FORWARD no iterations 04652 04653 //double Ftolerance = pow( d_macheps(), 0.5)*1.00; 04654 04655 04656 04657 //GZ out double Ftolerance = pow( d_macheps(), 0.5)*1000*KK; //Zhaohui UCD 10e6 for Pa, kg and m 1000 for kPa, ton and m 04658 04659 04660 04661 double Ftolerance = 1.0e-8; //Guanzhou gives an absolute tolerance 04662 04663 04664 04665 //opserr << Ftolerance << endlnn; 04666 04667 // double Ftolerance = pow(d_macheps(),(1.0/2.0))*1.0; 04668 04669 // double entry_kappa_cone = Criterion.kappa_cone_get(); 04670 04671 // double entry_kappa_cap = Criterion.kappa_cap_get(); 04672 04673 04674 04675 tensor aC(2, def_dim_2, 0.0); 04676 04677 stresstensor BEstress; 04678 04679 stresstensor residual; 04680 04681 tensor d2Qoverds2( 4, def_dim_4, 0.0); 04682 04683 tensor T( 4, def_dim_4, 0.0); 04684 04685 tensor Tinv( 4, def_dim_4, 0.0); 04686 04687 04688 04689 double Fold = 0.0; 04690 04691 tensor temp3lower; 04692 04693 tensor temp5; 04694 04695 double temp6 = 0.0; 04696 04697 double upper = 0.0; 04698 04699 04700 04701 stresstensor dsigma; 04702 04703 //stresstensor Dsigma; 04704 04705 stresstensor sigmaBack; 04706 04707 straintensor dPlasticStrain; // delta plastic strain 04708 04709 straintensor PlasticStrain; // Total plastic strain 04710 04711 straintensor incrPlasticStrain; 04712 04713 04714 04715 //double dq_ast = 0.0; // iterative change in internal variable (kappa in this case) 04716 04717 //double Dq_ast = 0.0; // incremental change in internal variable (kappa in this case) 04718 04719 //double q_ast = 0.0; // internal variable (kappa in this case) 04720 04721 //double q_ast_entry = 0.0; //internal variable from previous increment (kappa in this case) 04722 04723 04724 04725 int step_counter = 0; 04726 04727 //int MAX_STEP_COUNT = ; 04728 04729 // int MAX_STEP_COUNT = 0; 04730 04731 //int flag = 0; 04732 04733 04734 04735 straintensor strain_incr = strain_increment; 04736 04737 strain_incr.null_indices(); 04738 04739 stresstensor stress_increment = E("ijkl")*strain_incr("kl"); 04740 04741 stress_increment.null_indices(); 04742 04743 04744 04745 //stresstensor Return_stress; // stress to be returned can be either predictor 04746 04747 // or elastic plastic stress. 04748 04749 04750 04751 EPState ElasticPredictorEPS( startEPS ); 04752 04753 stresstensor elastic_predictor_stress = start_stress + stress_increment; 04754 04755 04756 04757 ElasticPredictorEPS.setStress( elastic_predictor_stress ); 04758 04759 // elastic_predictor_stress.reportshortpqtheta("\n . . . . ELASTIC PREDICTOR stress"); 04760 04761 04762 04763 //opserr.precision(5); 04764 04765 //opserr.width(10); 04766 04767 //opserr << "elastic predictor " << elastic_predictor_stress << endlnn; 04768 04769 04770 04771 stresstensor elastic_plastic_predictor_stress; 04772 04773 EPState EP_PredictorEPS( startEPS ); 04774 04775 04776 04777 //double f_start = material_point.YS->f( &startEPS ); 04778 04779 //opserr << " ************** f_start = " << f_start; 04780 04781 //::fprintf(stdout,"tst############## f_start = %.10e\n",f_start); 04782 04783 // f_pred = Criterion.f(elastic_predictor_stress); 04784 04785 //::fprintf(stdout,"tst############## f_pred = %.10e\n",f_pred); 04786 04787 //double f_start = getYS()->f( &startEPS ); 04788 04789 //opserr << "\n ************ Enter Backward f_star **** " << f_start; 04790 04791 04792 04793 double f_pred = getYS()->f( &ElasticPredictorEPS ); 04794 04795 //double f_pred = getYS()->f( &backwardEPS ); 04796 04797 //opserr << " **BE f_pred **** " << f_pred << endln; 04798 04799 //int region = 5; 04800 04801 04802 04803 //double h_s = 0.0; 04804 04805 //double xi_s = 0.0; 04806 04807 double hardMod_ = 0.0; 04808 04809 04810 04811 double h_s[4] = {0.0, 0.0, 0.0, 0.0}; 04812 04813 double xi_s[4] = {0.0, 0.0, 0.0, 0.0}; 04814 04815 stresstensor h_t[4]; 04816 04817 stresstensor xi_t[4]; 04818 04819 04820 04821 //region 04822 04823 //check for the region than proceede 04824 04825 //region = Criterion.influence_region(elastic_predictor_stress); 04826 04827 //if ( region == 1 ) // apex gray region 04828 04829 // { 04830 04831 // double pc_ = pc(); 04832 04833 // elastic_plastic_predictor_stress = 04834 04835 // elastic_plastic_predictor_stress.pqtheta2stress(pc_, 0.0, 0.0); 04836 04837 // return elastic_plastic_predictor_stress; 04838 04839 // } 04840 04841 04842 04843 if ( f_pred <= Ftolerance ) //Guanzhou changed 04844 04845 { 04846 04847 04848 04849 //Updating elastic strain increment 04850 04851 straintensor estrain = ElasticPredictorEPS.getElasticStrain(); //Guanzhou 04852 04853 straintensor tstrain = ElasticPredictorEPS.getStrain(); //Guanzhou 04854 04855 estrain = estrain + strain_increment; 04856 04857 tstrain = tstrain + strain_increment; 04858 04859 ElasticPredictorEPS.setElasticStrain( estrain ); 04860 04861 ElasticPredictorEPS.setStrain( tstrain ); 04862 04863 ElasticPredictorEPS.setdElasticStrain( strain_increment ); 04864 04865 //opserr<< "Elastic: Total strain" << tstrain << endln; 04866 04867 04868 04869 if ( getELT1() ) { 04870 04871 getELT1()->updateEeDm(&ElasticPredictorEPS, -st_vol, 0.0); 04872 04873 } 04874 04875 04876 04877 //Set Elasto-Plastic stiffness tensor 04878 04879 ElasticPredictorEPS.setEep(E); 04880 04881 ElasticPredictorEPS.setConverged(TRUE); 04882 04883 //E.printshort(" BE -- Eep "); 04884 04885 04886 04887 backwardEPS = ElasticPredictorEPS; 04888 04889 //opserr << "\n backwardEPS e " << backwardEPS.gete(); 04890 04891 04892 04893 //double p = (backwardEPS.getStress()).p_hydrostatic(); 04894 04895 //double ec = (backwardEPS.getec()) - (backwardEPS.getLam()) * log( p / (backwardEPS.getpo()) ); 04896 04897 //double st = (backwardEPS.getStrain()).Iinvariant1(); 04898 04899 //double pl_s = (backwardEPS.getPlasticStrain()).Iinvariant1(); 04900 04901 //double dpl_s = (backwardEPS.getdPlasticStrain()).Iinvariant1(); 04902 04903 //opserr << "E ec " << ec << " e " << backwardEPS.gete() << " psi " << (backwardEPS.gete() - ec) << " strain " << st << " t_pl " << pl_s << " d_pl " << dpl_s << "\n"; 04904 04905 return backwardEPS; 04906 04907 04908 04909 //opserr << "\nbackwardEPS" << backwardEPS; 04910 04911 //opserr << "\nElasticPredictorEPS " << ElasticPredictorEPS; 04912 04913 04914 04915 } 04916 04917 if ( f_pred > Ftolerance ) 04918 04919 { 04920 04921 04922 04923 //Let's put strict backwardEuler here 04924 04925 //Starting point by applying one Forward Euler step 04926 04927 EP_PredictorEPS = PredictorEPState( strain_incr); 04928 04929 //EP_PredictorEPS = ElasticPredictorEPS; 04930 04931 04932 04933 //opserr << " ----------Predictor Stress" << EP_PredictorEPS.getStress(); 04934 04935 //Setting the starting EPState with the starting internal vars in EPState 04936 04937 04938 04939 //MP->setEPS( EP_PredictorEPS ); 04940 04941 04942 04943 Felplpredictor = getYS()->f(&EP_PredictorEPS); 04944 04945 //opserr << " BE: F_elplpredictor " << Felplpredictor << endln; 04946 04947 04948 04949 04950 04951 //Kai absFelplpredictor = fabs(Felplpredictor); 04952 04953 if ( fabs(Felplpredictor) <= Ftolerance ) 04954 04955 { 04956 04957 //Forward Euler will do. 04958 04959 backwardEPS = EP_PredictorEPS; 04960 04961 return backwardEPS; //Guanzhou 04962 04963 //Return_stress = elastic_plastic_predictor_stress; 04964 04965 //Guanzhou out flag = 1; 04966 04967 } 04968 04969 04970 04971 //GZ out aC = getPS()->dQods( &EP_PredictorEPS ); 04972 04973 //GZ out dFods = getYS()->dFods( &EP_PredictorEPS ); 04974 04975 //GZ out dQods = getPS()->dQods( &EP_PredictorEPS ); 04976 04977 //GZ out 04978 04979 //GZ out temp2 = (dFods("ij")*E("ijkl"))*dQods("kl"); 04980 04981 //GZ out temp2.null_indices(); 04982 04983 //GZ out lower = temp2.trace(); 04984 04985 //GZ out 04986 04987 //GZ out //Guanzhou added internal evolution Mar2005 04988 04989 //GZ out hardMod_ = 0.0; 04990 04991 //GZ out tensor Hh(2, def_dim_2, 0.0); 04992 04993 //GZ out //Of 1st scalar internal vars 04994 04995 //GZ out if ( getELS1() ) { 04996 04997 //GZ out h_s[0] = getELS1()->h_s( &EP_PredictorEPS, getPS()); 04998 04999 //GZ out xi_s[0] = getYS()->xi_s1( &EP_PredictorEPS ); 05000 05001 //GZ out hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 05002 05003 //GZ out } 05004 05005 //GZ out 05006 05007 //GZ out //Of 2nd scalar internal vars 05008 05009 //GZ out if ( getELS2() ) { 05010 05011 //GZ out h_s[1] = getELS2()->h_s( &EP_PredictorEPS, getPS()); 05012 05013 //GZ out xi_s[1] = getYS()->xi_s2( &EP_PredictorEPS ); 05014 05015 //GZ out hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 05016 05017 //GZ out } 05018 05019 //GZ out 05020 05021 //GZ out //Of 3rd scalar internal vars 05022 05023 //GZ out if ( getELS3() ) { 05024 05025 //GZ out h_s[2] = getELS3()->h_s( &EP_PredictorEPS, getPS()); 05026 05027 //GZ out xi_s[2] = getYS()->xi_s3( &EP_PredictorEPS ); 05028 05029 //GZ out hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 05030 05031 //GZ out } 05032 05033 //GZ out 05034 05035 //GZ out //Of 4th scalar internal vars 05036 05037 //GZ out if ( getELS4() ) { 05038 05039 //GZ out h_s[3] = getELS4()->h_s( &EP_PredictorEPS, getPS()); 05040 05041 //GZ out xi_s[3] = getYS()->xi_s4( &EP_PredictorEPS ); 05042 05043 //GZ out hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 05044 05045 //GZ out } 05046 05047 //GZ out 05048 05049 //GZ out //Of tensorial internal var 05050 05051 //GZ out // 1st tensorial var 05052 05053 //GZ out if ( getELT1() ) { 05054 05055 //GZ out h_t[0] = getELT1()->h_t( &EP_PredictorEPS, getPS()); 05056 05057 //GZ out xi_t[0] = getYS()->xi_t1( &EP_PredictorEPS ); 05058 05059 //GZ out tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 05060 05061 //GZ out hardMod_ = hardMod_ + hm.trace(); 05062 05063 //GZ out } 05064 05065 //GZ out 05066 05067 //GZ out // 2nd tensorial var 05068 05069 //GZ out if ( getELT2() ) { 05070 05071 //GZ out h_t[1] = getELT2()->h_t( &EP_PredictorEPS, getPS()); 05072 05073 //GZ out xi_t[1] = getYS()->xi_t2( &EP_PredictorEPS ); 05074 05075 //GZ out tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 05076 05077 //GZ out hardMod_ = hardMod_ + hm.trace(); 05078 05079 //GZ out } 05080 05081 //GZ out 05082 05083 //GZ out // 3rd tensorial var 05084 05085 //GZ out if ( getELT3() ) { 05086 05087 //GZ out h_t[2] = getELT3()->h_t( &EP_PredictorEPS, getPS()); 05088 05089 //GZ out xi_t[2] = getYS()->xi_t3( &EP_PredictorEPS ); 05090 05091 //GZ out tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 05092 05093 //GZ out hardMod_ = hardMod_ + hm.trace(); 05094 05095 //GZ out } 05096 05097 //GZ out 05098 05099 //GZ out // 4th tensorial var 05100 05101 //GZ out if ( getELT4() ) { 05102 05103 //GZ out h_t[3] = getELT4()->h_t( &EP_PredictorEPS, getPS()); 05104 05105 //GZ out xi_t[3] = getYS()->xi_t4( &EP_PredictorEPS ); 05106 05107 //GZ out tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 05108 05109 //GZ out hardMod_ = hardMod_ + hm.trace(); 05110 05111 //GZ out } 05112 05113 //GZ out 05114 05115 //GZ out 05116 05117 //GZ out lower = lower - hardMod_; 05118 05119 //GZ out 05120 05121 //GZ out 05122 05123 //GZ out Delta_lambda = f_pred/lower; //Guanzhou Mar2005 05124 05125 //GZ out //::printf(" Delta_lambda = f_pred/lower = %.8e\n", Delta_lambda); 05126 05127 //GZ out //// Delta_lambda = Felplpredictor/lower; 05128 05129 //GZ out ////::printf(" Delta_lambda = Felplpredictor/lower =%.8e \n", Delta_lambda); 05130 05131 //GZ out 05132 05133 //GZ out elastic_plastic_predictor_stress = elastic_predictor_stress - E("ijkl")*aC("kl")*Delta_lambda; 05134 05135 //GZ out elastic_plastic_predictor_stress.null_indices(); 05136 05137 //GZ out EP_PredictorEPS.setStress( elastic_plastic_predictor_stress ); 05138 05139 //GZ out 05140 05141 //GZ out //Guanzhou, update internal 05142 05143 //GZ out 05144 05145 //GZ out int NS = EP_PredictorEPS.getNScalarVar(); 05146 05147 //GZ out int NT = EP_PredictorEPS.getNTensorVar(); 05148 05149 //GZ out 05150 05151 //GZ out double dS = 0; 05152 05153 //GZ out double S = 0; 05154 05155 //GZ out //double new_S = 0; 05156 05157 //GZ out 05158 05159 //GZ out stresstensor dT; 05160 05161 //GZ out stresstensor Tv; 05162 05163 //GZ out stresstensor new_T; 05164 05165 //GZ out 05166 05167 //GZ out int ii; 05168 05169 //GZ out for (ii = 1; ii <= NS; ii++) { 05170 05171 //GZ out dS = Delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var 05172 05173 //GZ out S = EP_PredictorEPS.getScalarVar(ii); // Get the old value of the scalar internal var 05174 05175 //GZ out EP_PredictorEPS.setScalarVar(ii, S + dS ); // Update internal scalar var 05176 05177 //GZ out } 05178 05179 //GZ out 05180 05181 //GZ out for (ii = 1; ii <= NT; ii++) { 05182 05183 //GZ out dT = h_t[ii-1] * Delta_lambda; // Increment to the tensor internal var 05184 05185 //GZ out Tv = EP_PredictorEPS.getTensorVar(ii); // Get the old value of the tensor internal var 05186 05187 //GZ out new_T = Tv + dT; 05188 05189 //GZ out EP_PredictorEPS.setTensorVar(ii, new_T ); // Update tensorial scalar var 05190 05191 //GZ out } 05192 05193 //GZ out 05194 05195 //GZ out Felplpredictor = getYS()->f(&EP_PredictorEPS); 05196 05197 05198 05199 //Guanzhou out Mar2005 //Zhaohui modified, sometimes give much better convergence rate 05200 05201 //Guanzhou out Mar2005 elastic_plastic_predictor_stress = EP_PredictorEPS.getStress(); 05202 05203 //Guanzhou out Mar2005 //opserr << "elastic_plastic_predictor_stress" << elastic_plastic_predictor_stress; 05204 05205 05206 05207 //opserr.precision(5); 05208 05209 //opserr.width(10); 05210 05211 //opserr << " " << EP_PredictorEPS.getStress().p_hydrostatic() << " "; 05212 05213 05214 05215 //opserr.precision(5); 05216 05217 //opserr.width(10); 05218 05219 //opserr << EP_PredictorEPS.getStress().q_deviatoric()<< " "; 05220 05221 05222 05223 //opserr.precision(5); 05224 05225 //opserr.width(10); 05226 05227 //opserr << Delta_lambda << endlnn; 05228 05229 05230 05231 //elastic_plastic_predictor_stress.reportshort("......elastic_plastic_predictor_stress"); 05232 05233 //::printf(" F(elastic_plastic_predictor_stress) = %.8e\n", Criterion.f(elastic_plastic_predictor_stress)); 05234 05235 05236 05237 //h_s = MP->ELS1->h_s( &EP_PredictorEPS, MP->PS ); 05238 05240 05242 05243 05244 05245 //q_ast_entry = Criterion.kappa_get(elastic_plastic_predictor_stress); 05246 05247 //Dq_ast = Delta_lambda * h_ * just_this_PP; 05248 05249 //q_ast = q_ast_entry + Dq_ast; 05250 05251 05252 05253 //Zhaohui comments: internal vars are alreary evolued in ForwardEulerEPS(...), not necessary here!!! 05254 05255 //..dS = Delta_lambda * h_s ; // Increment to the internal var 05256 05257 //..S = EP_PredictorEPS.getScalarVar(1); // Get the old value of the internal var 05258 05259 //..new_S = S + dS; 05260 05261 //..opserr << "Internal Var : " << new_S << endlnn; 05262 05263 //..EP_PredictorEPS.setScalarVar(1, new_S); // Get the old value of the internal var 05264 05265 05266 05267 //::fprintf(stdout," ######## predictor --> Dq_ast=%.8e q_ast=%.8e\n", Dq_ast, q_ast); 05268 05269 //::fprintf(stderr," ######## predictor --> Dq_ast=%.8e q_ast=%.8e\n", Dq_ast, q_ast); 05270 05271 05272 05273 //Criterion.kappa_set( sigmaBack, q_ast); //???? 05274 05275 //current_lambda_set(Delta_lambda); //???? 05276 05277 05278 05279 //::printf(" Delta_lambda = %.8e\n", Delta_lambda); 05280 05281 //::printf("step = pre iteracija #############################-- q_ast = %.10e \n", q_ast); 05282 05283 //::printf("step = pre iteracija posle predictora ###########-- Dq_ast = %.10e \n",Dq_ast); 05284 05285 //********** 05286 05287 //********** 05288 05289 //::printf("\nDelta_lambda before BE = %.10e \n", Delta_lambda ); 05290 05291 05292 05293 //========================== main part of iteration ======================= 05294 05295 // while ( absFelplpredictor > Ftolerance && 05296 05297 double Felplold = 0.0; 05298 05299 Delta_lambda = EP_PredictorEPS.Delta_lambda; 05300 05301 while (( fabs(fabs(Felplpredictor) - fabs(Felplold)) > 1.0e-6 ) 05302 05303 && ( fabs(fabs(Felplpredictor) ) > 1.0e-5 ) 05304 05305 && (step_counter < MAX_STEP_COUNT) )// if more iterations than prescribed 05306 05307 //while ( ( fabs(fabs(Felplpredictor) ) > Ftolerance ) && (step_counter < MAX_STEP_COUNT) )// if more iterations than prescribed 05308 05309 //out07may97 do 05310 05311 { 05312 05313 //opserr << "Iteration " << step_counter << " F " << Felplpredictor << " delta " << Delta_lambda << endln; 05314 05315 aC = getPS()->dQods( &EP_PredictorEPS ); 05316 05317 BEstress = elastic_predictor_stress - E("ijkl")*aC("kl")*Delta_lambda; 05318 05319 //BEstress.reportshort("......BEstress "); 05320 05322 05323 BEstress.null_indices(); 05324 05325 // Felplpredictor = Criterion.f(BEstress); 05326 05327 // ::printf("\nF_backward_Euler BE = %.10e \n", Felplpredictor); 05328 05329 05330 05331 //Guanzhou out Mar2005 residual = elastic_plastic_predictor_stress - BEstress; 05332 05333 residual = EP_PredictorEPS.getStress() - BEstress; 05334 05335 05336 05337 //residual.reportshortpqtheta("\n......residual "); 05338 05339 // double ComplementaryEnergy = (residual("ij")*D("ijkl")*residual("ij")).trace(); 05340 05341 //::printf("\n Residual ComplementaryEnergy = %.16e\n", ComplementaryEnergy); 05342 05343 05344 05346 05347 05348 05349 05350 05351 //d2Qoverds2 = Criterion.d2Qods2(elastic_plastic_predictor_stress); 05352 05353 d2Qoverds2 = getPS()->d2Qods2( &EP_PredictorEPS ); 05354 05355 //d2Qoverds2.print(); 05356 05357 05358 05359 T = I_ikjl + E("ijkl")*d2Qoverds2("klmn")*Delta_lambda; 05360 05361 T.null_indices(); 05362 05363 05364 05365 Tinv = T.inverse(); 05366 05367 05368 05369 //Guanzhou out Mar2005 05370 05371 //Guanzhou out Mar2005 //Z Cheng add: H 05372 05373 //Guanzhou out Mar2005 tensor H( 2, def_dim_2, 0.0); 05374 05375 //Guanzhou out Mar2005 H = dQods; 05376 05377 05378 05379 //dFods = Criterion.dFods(elastic_plastic_predictor_stress); 05380 05381 //dQods = Criterion.dQods(elastic_plastic_predictor_stress); 05382 05383 dFods = getYS()->dFods( &EP_PredictorEPS ); 05384 05385 dQods = getPS()->dQods( &EP_PredictorEPS ); 05386 05387 05388 05389 //Guanzhou Mar2005 05390 05391 tensor H( 2, def_dim_2, 0.0); 05392 05393 H = dQods; 05394 05395 05396 05397 //Fold = Criterion.f(elastic_plastic_predictor_stress); 05398 05399 Fold = getYS()->f( &EP_PredictorEPS ); 05400 05401 05402 05403 lower = 0.0; // this is old temp variable used here again :-) 05404 05405 //h_ = h(elastic_plastic_predictor_stress); 05406 05407 //xi_ = xi(elastic_plastic_predictor_stress); 05408 05409 05410 05411 //h_s = MP->ELS1->h_s( &EP_PredictorEPS, MP->PS ); 05412 05413 //xi_s = MP->YS->xi_s1( &EP_PredictorEPS ); 05414 05415 //hardMod_ = h_s * xi_s; 05416 05417 05418 05419 // Evaluating the hardening modulus: sum of (df/dq*) * qbar 05420 05421 05422 05423 hardMod_ = 0.0; 05424 05425 tensor Hh(2, def_dim_2, 0.0); 05426 05427 //Of 1st scalar internal vars 05428 05429 05430 05431 if ( getELS1() ) { 05432 05433 h_s[0] = getELS1()->h_s( &EP_PredictorEPS, getPS()); 05434 05435 xi_s[0] = getYS()->xi_s1( &EP_PredictorEPS ); 05436 05437 hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 05438 05439 Hh = Hh + getPS()->d2Qodsds1( &EP_PredictorEPS ) *h_s[0] *Delta_lambda; 05440 05441 } 05442 05443 05444 05445 //Of 2nd scalar internal vars 05446 05447 if ( getELS2() ) { 05448 05449 h_s[1] = getELS2()->h_s( &EP_PredictorEPS, getPS()); 05450 05451 xi_s[1] = getYS()->xi_s2( &EP_PredictorEPS ); 05452 05453 hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 05454 05455 Hh = Hh + getPS()->d2Qodsds2( &EP_PredictorEPS ) *h_s[1] *Delta_lambda; 05456 05457 } 05458 05459 05460 05461 //Of 3rd scalar internal vars 05462 05463 if ( getELS3() ) { 05464 05465 h_s[2] = getELS3()->h_s( &EP_PredictorEPS, getPS()); 05466 05467 xi_s[2] = getYS()->xi_s3( &EP_PredictorEPS ); 05468 05469 hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 05470 05471 Hh = Hh + getPS()->d2Qodsds3( &EP_PredictorEPS ) *h_s[2] *Delta_lambda; 05472 05473 } 05474 05475 05476 05477 //Of 4th scalar internal vars 05478 05479 if ( getELS4() ) { 05480 05481 h_s[3] = getELS4()->h_s( &EP_PredictorEPS, getPS()); 05482 05483 xi_s[3] = getYS()->xi_s4( &EP_PredictorEPS ); 05484 05485 hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 05486 05487 Hh = Hh + getPS()->d2Qodsds4( &EP_PredictorEPS ) *h_s[3] *Delta_lambda; 05488 05489 } 05490 05491 05492 05493 //Of tensorial internal var 05494 05495 // 1st tensorial var 05496 05497 if ( getELT1() ) { 05498 05499 h_t[0] = getELT1()->h_t( &EP_PredictorEPS, getPS()); 05500 05501 xi_t[0] = getYS()->xi_t1( &EP_PredictorEPS ); 05502 05503 tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 05504 05505 hardMod_ = hardMod_ + hm.trace(); 05506 05507 Hh = Hh + (getPS()->d2Qodsdt1( &EP_PredictorEPS ))("ijmn") *(h_t[0])("mn") *Delta_lambda; 05508 05509 } 05510 05511 05512 05513 // 2nd tensorial var 05514 05515 if ( getELT2() ) { 05516 05517 h_t[1] = getELT2()->h_t( &EP_PredictorEPS, getPS()); 05518 05519 xi_t[1] = getYS()->xi_t2( &EP_PredictorEPS ); 05520 05521 tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 05522 05523 hardMod_ = hardMod_ + hm.trace(); 05524 05525 Hh = Hh + (getPS()->d2Qodsdt2( &EP_PredictorEPS ))("ijmn") *(h_t[1])("mn") *Delta_lambda; 05526 05527 } 05528 05529 05530 05531 // 3rd tensorial var 05532 05533 if ( getELT3() ) { 05534 05535 h_t[2] = getELT3()->h_t( &EP_PredictorEPS, getPS()); 05536 05537 xi_t[2] = getYS()->xi_t3( &EP_PredictorEPS ); 05538 05539 tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 05540 05541 hardMod_ = hardMod_ + hm.trace(); 05542 05543 Hh = Hh + (getPS()->d2Qodsdt3( &EP_PredictorEPS ))("ijmn") *(h_t[2])("mn") *Delta_lambda; 05544 05545 } 05546 05547 05548 05549 // 4th tensorial var 05550 05551 if ( getELT4() ) { 05552 05553 h_t[3] = getELT4()->h_t( &EP_PredictorEPS, getPS()); 05554 05555 xi_t[3] = getYS()->xi_t4( &EP_PredictorEPS ); 05556 05557 tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 05558 05559 hardMod_ = hardMod_ + hm.trace(); 05560 05561 Hh = Hh + (getPS()->d2Qodsdt4( &EP_PredictorEPS ))("ijmn") *(h_t[3])("mn") *Delta_lambda; 05562 05563 } 05564 05565 05566 05567 // Subtract accumulated hardMod_ from lower 05568 05569 //lower = lower - hardMod_; 05570 05571 05572 05573 //hardMod_ = hardMod_ * just_this_PP; 05574 05575 //::printf("\n BackwardEulerStress .. hardMod_ = %.10e \n", hardMod_ ); 05576 05577 //outfornow d2Qodqast = d2Qoverdqast(elastic_plastic_predictor_stress); 05578 05579 //outfornow dQodsextended = dQods + d2Qodqast * Delta_lambda * h_; 05580 05581 //outfornow temp3lower = dFods("mn")*Tinv("ijmn")*E("ijkl")*dQodsextended("kl"); 05582 05583 // temp3lower = dFods("mn")*Tinv("ijmn")*E("ijkl")*dQods("kl"); // Bug found, Z Cheng, Jan 2004 05584 05585 Hh.null_indices(); 05586 05587 05588 05589 H = H + Hh; 05590 05591 05592 05593 temp3lower = dFods("mn")*Tinv("ijmn")*E("ijkl")*H("kl"); 05594 05595 temp3lower.null_indices(); 05596 05597 lower = temp3lower.trace(); 05598 05599 lower = lower - hardMod_; 05600 05601 //opserr << " 1st hardMod_ " << hardMod_ << "\n"; 05602 05603 05604 05605 temp5 = (residual("ij") * Tinv("ijmn")) * dFods("mn"); 05606 05607 temp6 = temp5.trace(); 05608 05609 //The same as the above but more computation 05610 05611 //temp5 = dFods("mn") * residual("ij") * Tinv("ijmn"); 05612 05613 //temp6 = temp5.trace(); 05614 05615 05616 05617 upper = Fold - temp6; 05618 05619 05620 05621 //================================================================================ 05622 05623 //dlambda 05624 05625 delta_lambda = upper / lower; 05626 05627 //opserr << "upper " << upper << " lower " << lower << endln; 05628 05629 Delta_lambda = Delta_lambda + delta_lambda; 05630 05631 if ( Delta_lambda < 0.0 ) Delta_lambda = 0.0;//Guanzhou 05632 05633 05634 05635 //Guanzhou corrected May2005 // Zhaohui_____10-01-2000 not sure xxxxxxxxxxxxxxxxxxx 05636 05637 dPlasticStrain = dQods("kl") * delta_lambda; 05638 05639 incrPlasticStrain = incrPlasticStrain + dPlasticStrain; 05640 05641 //Updating elastic strain increment 05642 05643 straintensor estrain = EP_PredictorEPS.getElasticStrain(); //Guanzhou 05644 05645 straintensor tstrain = EP_PredictorEPS.getStrain(); //Guanzhou 05646 05647 straintensor pstrain = tstrain - estrain; 05648 05649 pstrain = pstrain + incrPlasticStrain; 05650 05651 estrain = tstrain - pstrain; 05652 05653 EP_PredictorEPS.setElasticStrain( estrain ); 05654 05655 EP_PredictorEPS.setPlasticStrain( pstrain ); 05656 05657 EP_PredictorEPS.setdPlasticStrain( incrPlasticStrain ); 05658 05659 straintensor incrElasticStrain = incrPlasticStrain * (-1.0); 05660 05661 EP_PredictorEPS.setdElasticStrain( incrElasticStrain ); 05662 05663 EP_PredictorEPS.Delta_lambda = Delta_lambda; 05664 05665 05666 05667 //opserr<< "Elastic: Total strain" << tstrain << endln; 05668 05669 05670 05671 int err; 05672 05673 if ( getELT1() ) { 05674 05675 double pl_st_vol = pstrain.Iinvariant1(); //Joey 02-17-03 05676 05677 // D > 0 compressive -> Iinv > 0 -> de < 0 correct! 05678 05679 err = getELT1()->updateEeDm(&EP_PredictorEPS, pl_st_vol, Delta_lambda); 05680 05681 } 05682 05683 05684 05685 //::printf(" >> %d Delta_lambda = %.8e", step_counter, Delta_lambda); 05686 05687 // stari umesto dQodsextended za stari = dQods 05688 05689 //outfornow dsigma = 05690 05691 //outfornow ((residual("ij")*Tinv("ijmn"))+ 05692 05693 //outfornow ((E("ijkl")*dQodsextended("kl"))*Tinv("ijmn")*Delta_lambda) )*-1.0; 05694 05695 //::printf(" Delta_lambda = %.8e\n", Delta_lambda); 05696 05697 //dsigma = ( (residual("ij")*Tinv("ijmn") )+ 05698 05699 // ( (E("ijkl")*dQods("kl"))*Tinv("ijmn")*delta_lambda) )*(-1.0); //*-1.0??? // Bug found, Z Cheng, Jan 2004 05700 05701 dsigma = ( (residual("ij")*Tinv("ijmn") )+ 05702 05703 ( (E("ijkl")*H("kl"))*Tinv("ijmn")*delta_lambda) )*(-1.0); //Guanzhou fixed Apr2005 05704 05705 dsigma.null_indices(); 05706 05707 05708 05709 sigmaBack = EP_PredictorEPS.getStress() + dsigma; 05710 05711 EP_PredictorEPS.setStress(sigmaBack); 05712 05713 05714 05715 //dsigma.reportshortpqtheta("\n......dsigma "); 05716 05717 //::printf(" ......... in NR loop delta_lambda = %.16e\n", delta_lambda); 05718 05719 //::printf(" ......... in NR loop Delta_lambda = %.16e\n", Delta_lambda); 05720 05721 05722 05723 //dq_ast = delta_lambda * h_ * just_this_PP; 05724 05725 //Dq_ast += dq_ast; 05726 05727 05728 05729 //dS = delta_lambda * h_s ; // Increment to the internal var 05730 05731 //S = EP_PredictorEPS.getScalarVar(1); // Get the old value of the internal var 05732 05733 //new_S = S + dS; 05734 05735 //EP_PredictorEPS.setScalarVar(1, new_S); 05736 05737 05738 05739 //Evolve the surfaces and hardening vars 05740 05741 int NS = EP_PredictorEPS.getNScalarVar(); 05742 05743 int NT = EP_PredictorEPS.getNTensorVar(); 05744 05745 05746 05747 double dS = 0; 05748 05749 double S = 0; 05750 05751 //double new_S = 0; 05752 05753 05754 05755 stresstensor dT; 05756 05757 stresstensor T; 05758 05759 stresstensor new_T; 05760 05761 05762 05763 //if ( delta_lambda < 0) delta_lambda = 0; 05764 05765 int ii; 05766 05767 for (ii = 1; ii <= NS; ii++) { 05768 05769 dS = delta_lambda * h_s[ii-1] ; //Guanzhou fixed Apr2005 // Increment to the scalar internal var 05770 05771 S = EP_PredictorEPS.getScalarVar(ii); //Guanzhou fixed Apr2005// Get the old value of the scalar internal var 05772 05773 EP_PredictorEPS.setScalarVar(ii, S + dS ); // Update internal scalar var 05774 05775 } 05776 05777 05778 05779 for (ii = 1; ii <= NT; ii++) { 05780 05781 dT = h_t[ii-1] * delta_lambda; //Guanzhou fixed Apr2005 // Increment to the tensor internal var 05782 05783 T = EP_PredictorEPS.getTensorVar(ii); //Guanzhou fixed Apr2005// Get the old value of the tensor internal var 05784 05785 new_T = T + dT; 05786 05787 EP_PredictorEPS.setTensorVar(ii, new_T ); // Update tensorial scalar var 05788 05789 } 05790 05791 05792 05793 Felplold = Felplpredictor; 05794 05795 Felplpredictor = getYS()->f( &EP_PredictorEPS ); 05796 05797 05798 05799 //======= Dq_ast = Delta_lambda * h_ * just_this_PP; 05800 05801 //q_ast = q_ast_entry + Dq_ast; 05802 05803 05804 05805 //::fprintf(stdout," ######## step = %3d --> Dq_ast=%.8e q_ast=%.8e\n", 05806 05807 // step_counter, Dq_ast, q_ast); 05808 05809 //::fprintf(stderr," ######## step = %3d --> Dq_ast=%.8e q_ast=%.8e\n", 05810 05811 // step_counter, Dq_ast, q_ast); 05812 05813 05814 05815 //current_lambda_set(Delta_lambda); 05816 05817 //.... elastic_plastic_predictor_stress.reportshort("elplpredstress"); 05818 05819 //.... dsigma.reportshort("dsigma"); 05820 05821 05822 05823 //sigmaBack.reportshortpqtheta("\n before======== SigmaBack"); 05824 05825 //Guanzhou out Mar2005 sigmaBack = elastic_plastic_predictor_stress + dsigma; 05826 05827 //sigmaBack.deviator().reportshort("\n after ======== SigmaBack"); 05828 05829 //sigmaBack.reportshortpqtheta("\n after ======== SigmaBack"); 05830 05831 05832 05833 //Guanzhou out Mar2005 //temp trick 05834 05835 //Guanzhou out Mar2005 if (sigmaBack.p_hydrostatic() > 0) 05836 05837 //Guanzhou out Mar2005 { 05838 05839 //Guanzhou out Mar2005 //====== sigmaBack = elastic_predictor_stress + Dsigma; 05840 05841 //Guanzhou out Mar2005 //sigmaBack.reportshortpqtheta("BE................ NR sigmaBack "); 05842 05843 //Guanzhou out Mar2005 //sigmaBack.reportAnim(); 05844 05845 //Guanzhou out Mar2005 //::fprintf(stdout,"Anim BEpoint0%d = {Sin[theta]*q, p, Cos[theta]*q} \n",step_counter+1); 05846 05847 //Guanzhou out Mar2005 ////::fprintf(stdout,"Anim BEpoint0%dP = Point[BEpoint0%d] \n",step_counter+1,step_counter+1); 05848 05849 //Guanzhou out Mar2005 //::fprintf(stdout,"Anim \n"); 05850 05851 //Guanzhou out Mar2005 05852 05853 //Guanzhou out Mar2005 //Criterion.kappa_set( sigmaBack, q_ast) ; 05854 05855 //Guanzhou out Mar2005 EP_PredictorEPS.setStress( sigmaBack ); 05856 05857 //Guanzhou out Mar2005 05858 05859 //Guanzhou out Mar2005 //Felplpredictor = Criterion.f(sigmaBack); 05860 05861 //Guanzhou out Mar2005 //Kai absFelplpredictor = fabs(Felplpredictor); 05862 05863 //Guanzhou out Mar2005 //::printf(" F_bE=%.10e (%.10e)\n", Felplpredictor,Ftolerance); 05864 05865 //Guanzhou out Mar2005 Felplpredictor = getYS()->f( &EP_PredictorEPS ); 05866 05867 //Guanzhou out Mar2005 //::printf(" F_BE: step=%5d F= %.10e (%.10e)\n", step_counter, Felplpredictor, Ftolerance); 05868 05869 //Guanzhou out Mar2005 } 05870 05871 //Guanzhou out Mar2005 else 05872 05873 //Guanzhou out Mar2005 { 05874 05875 //Guanzhou out Mar2005 sigmaBack= sigmaBack.pqtheta2stress(0.1, 0.0, 0.0); 05876 05877 //Guanzhou out Mar2005 Felplpredictor = 0; 05878 05879 //Guanzhou out Mar2005 backwardEPS.setStress(sigmaBack); 05880 05881 //Guanzhou out Mar2005 backwardEPS.setConverged(TRUE); 05882 05883 //Guanzhou out Mar2005 return backwardEPS; 05884 05885 //Guanzhou out Mar2005 } 05886 05887 05888 05889 //double tempkappa1 = kappa_cone_get(); 05890 05891 //double tempdFodeta = dFoverdeta(sigmaBack); 05892 05893 //double tempdetaodkappa = detaoverdkappa(tempkappa1); 05894 05895 //::printf(" h_=%.6e xi_=%.6e, dFodeta=%.6e, detaodkappa=%.6e, hardMod_=%.6e\n", 05896 05897 // h_, xi_,tempdFodeta, tempdetaodkappa, hardMod_); 05898 05899 //::printf(" upper = %.6e lower = %.6e\n", upper, lower); 05900 05901 //::printf(" q_ast_entry=%.6e Dq_ast=%.6e Delta_lambda=%.6e\n", 05902 05903 // q_ast_entry, Dq_ast, Delta_lambda); 05904 05905 05906 05907 // now prepare new step 05908 05909 //Guanzhou out Mar2005 elastic_plastic_predictor_stress = sigmaBack; 05910 05911 05912 05913 //Output for plotting 05914 05915 //opserr.precision(5); 05916 05917 //opserr.width(10); 05918 05919 //opserr << " " << sigmaBack.p_hydrostatic() << " "; 05920 05921 //opserr << " " << sigmaBack.p_hydrostatic() << " "; 05922 05923 05924 05925 //opserr.precision(5); 05926 05927 //opserr.width(10); 05928 05929 //opserr << sigmaBack.q_deviatoric() << " "; 05930 05931 //opserr << sigmaBack.q_deviatoric() << " "; 05932 05933 05934 05935 //opserr.precision(5); 05936 05937 //opserr.width(10); 05938 05939 //opserr << " Felpl " << Felplpredictor; 05940 05941 //opserr << " Felpl " << Felplpredictor; 05942 05943 05944 05945 //opserr.precision(5); 05946 05947 //opserr.width(10); 05948 05949 //opserr << " tol " << Ftolerance << endlnn; 05950 05951 //opserr << " tol " << Ftolerance << " " << step_counter << endlnn; 05952 05953 05954 05955 //::printf(" ........................... end of step %d\n", step_counter);// getchar(); 05956 05957 step_counter++; 05958 05959 05960 05961 //int err = 0; 05962 05963 //if ( getELT1() ) { 05964 05965 // double pl_st_vol = dPlasticStrain.Iinvariant1(); //Joey 02-17-03 Iinvariant1 > 0 for compression 05966 05967 // opserr << " Pl update " << pl_st_vol << "\n"; 05968 05969 // err = getELT1()->updateEeDm(&EP_PredictorEPS, - pl_st_vol, delta_lambda); // D > 0 --> contraction 05970 05971 //} 05972 05973 05974 05975 } 05976 05977 05978 05979 //opserr << " " << sigmaBack.p_hydrostatic() << " "; 05980 05981 //opserr << sigmaBack.q_deviatoric() << " "; 05982 05983 //opserr << " Felpl " << Felplpredictor; 05984 05985 //opserr << " tol " << Ftolerance << " " << step_counter << endln; 05986 05987 05988 05990 05991 if ( step_counter >= MAX_STEP_COUNT ) 05992 05993 { 05994 05995 //g3ErrorHandler->warning("Template3Dep::BackwardEulerEPState Step_counter > MAX_STEP_COUNT %d iterations", MAX_STEP_COUNT ); 05996 05997 //Guanzhou 05998 05999 opserr << "Template3Dep::BackwardEuler(), failed to converge in " << step_counter << "steps!!!" << '\n'; 06000 06001 EP_PredictorEPS.setConverged( false ); 06002 06003 exit(1); 06004 06005 //::exit(1); 06006 06007 } 06008 06009 06010 06011 //Guanzhou out int err = 0; 06012 06013 //Guanzhou out if ( getELT1() ) { 06014 06015 //Guanzhou out double pl_st_vol = PlasticStrain.Iinvariant1(); //Joey 02-17-03 06016 06017 //Guanzhou out // D > 0 compressive -> Iinv > 0 -> de < 0 correct! 06018 06019 //Guanzhou out err = getELT1()->updateEeDm(&EP_PredictorEPS, pl_st_vol, Delta_lambda); 06020 06021 //Guanzhou out } 06022 06023 06024 06025 //out07may97 while ( absFelplpredictor > Ftolerance && 06026 06027 //out07may97 step_counter <= MAX_STEP_COUNT ); // if more iterations than prescribed 06028 06029 06030 06031 //********** 06032 06033 //********** 06034 06035 //********** 06036 06037 //********** 06038 06039 06040 06041 // already set everything 06042 06043 // Need to genarate Eep and set strains and stresses 06044 06045 //if ( ( flag !=1) && (step_counter < MAX_STEP_COUNT) ) 06046 06047 //Guanzhou out Mar2005 if ( ( flag !=1) ) { 06048 06049 //Guanzhou out Mar2005 06050 06051 06052 06053 //Return_stress = elastic_plastic_predictor_stress; 06054 06055 //Criterion.kappa_set( Return_stress, q_ast) ; 06056 06057 06058 06059 // Generating Consistent Stiffness Tensor Eep 06060 06061 tensor I2("I", 2, def_dim_2); 06062 06063 tensor I_ijkl = I2("ij")*I2("kl"); 06064 06065 I_ijkl.null_indices(); 06066 06067 tensor I_ikjl = I_ijkl.transpose0110(); 06068 06069 06070 06071 06072 06073 dQods = getPS()->dQods( &EP_PredictorEPS ); // this is m_ij 06074 06075 tensor temp2 = E("ijkl")*dQods("kl"); 06076 06077 temp2.null_indices(); 06078 06079 dFods = getYS()->dFods( &EP_PredictorEPS ); // this is n_ij 06080 06081 d2Qoverds2 = getPS()->d2Qods2( &EP_PredictorEPS ); 06082 06083 06084 06085 tensor T = I_ikjl + E("ijkl")*d2Qoverds2("klmn")*Delta_lambda; 06086 06087 //tensor tt = E("ijkl")*d2Qoverds2("klmn")*Delta_lambda; 06088 06089 //tt.printshort("temp tt"); 06090 06091 //T = I_ikjl + tt; 06092 06093 T.null_indices(); 06094 06095 tensor Tinv = T.inverse(); 06096 06097 06098 06099 tensor R = Tinv("ijmn")*E("ijkl"); 06100 06101 R.null_indices(); 06102 06103 06104 06105 tensor Hh(2, def_dim_2, 0.0); 06106 06107 //Hh.Reset_to(0.0); 06108 06109 tensor H(2, def_dim_2, 0.0); 06110 06111 H =dQods; 06112 06113 06114 06115 //Of 1st scalar internal vars 06116 06117 if ( getELS1() ) { 06118 06119 h_s[0] = getELS1()->h_s( &EP_PredictorEPS, getPS()); 06120 06121 xi_s[0] = getYS()->xi_s1( &EP_PredictorEPS ); 06122 06123 hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 06124 06125 Hh = Hh + getPS()->d2Qodsds1( &EP_PredictorEPS ) *h_s[0] *Delta_lambda; 06126 06127 } 06128 06129 06130 06131 //Of 2nd scalar internal vars 06132 06133 if ( getELS2() ) { 06134 06135 h_s[1] = getELS2()->h_s( &EP_PredictorEPS, getPS()); 06136 06137 xi_s[1] = getYS()->xi_s2( &EP_PredictorEPS ); 06138 06139 hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 06140 06141 Hh = Hh + getPS()->d2Qodsds2( &EP_PredictorEPS ) *h_s[1] *Delta_lambda; 06142 06143 } 06144 06145 06146 06147 //Of 3rd scalar internal vars 06148 06149 if ( getELS3() ) { 06150 06151 h_s[2] = getELS3()->h_s( &EP_PredictorEPS, getPS()); 06152 06153 xi_s[2] = getYS()->xi_s3( &EP_PredictorEPS ); 06154 06155 hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 06156 06157 Hh = Hh + getPS()->d2Qodsds3( &EP_PredictorEPS ) *h_s[2] *Delta_lambda; 06158 06159 } 06160 06161 06162 06163 //Of 4th scalar internal vars 06164 06165 if ( getELS4() ) { 06166 06167 h_s[3] = getELS4()->h_s( &EP_PredictorEPS, getPS()); 06168 06169 xi_s[3] = getYS()->xi_s4( &EP_PredictorEPS ); 06170 06171 hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 06172 06173 Hh = Hh + getPS()->d2Qodsds4( &EP_PredictorEPS ) *h_s[3] *Delta_lambda; 06174 06175 } 06176 06177 06178 06179 //Of tensorial internal var 06180 06181 // 1st tensorial var 06182 06183 if ( getELT1() ) { 06184 06185 h_t[0] = getELT1()->h_t( &EP_PredictorEPS, getPS()); 06186 06187 xi_t[0] = getYS()->xi_t1( &EP_PredictorEPS ); 06188 06189 tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 06190 06191 hardMod_ = hardMod_ + hm.trace(); 06192 06193 Hh = Hh + (getPS()->d2Qodsdt1( &EP_PredictorEPS ))("ijmn") *(h_t[0])("mn") *Delta_lambda; 06194 06195 } 06196 06197 06198 06199 // 2nd tensorial var 06200 06201 if ( getELT2() ) { 06202 06203 h_t[1] = getELT2()->h_t( &EP_PredictorEPS, getPS()); 06204 06205 xi_t[1] = getYS()->xi_t2( &EP_PredictorEPS ); 06206 06207 tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 06208 06209 hardMod_ = hardMod_ + hm.trace(); 06210 06211 Hh = Hh + (getPS()->d2Qodsdt2( &EP_PredictorEPS ))("ijmn") *(h_t[1])("mn") *Delta_lambda; 06212 06213 } 06214 06215 06216 06217 // 3rd tensorial var 06218 06219 if ( getELT3() ) { 06220 06221 h_t[2] = getELT3()->h_t( &EP_PredictorEPS, getPS()); 06222 06223 xi_t[2] = getYS()->xi_t3( &EP_PredictorEPS ); 06224 06225 tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 06226 06227 hardMod_ = hardMod_ + hm.trace(); 06228 06229 Hh = Hh + (getPS()->d2Qodsdt3( &EP_PredictorEPS ))("ijmn") *(h_t[2])("mn") *Delta_lambda; 06230 06231 } 06232 06233 06234 06235 // 4th tensorial var 06236 06237 if ( getELT4() ) { 06238 06239 h_t[3] = getELT4()->h_t( &EP_PredictorEPS, getPS()); 06240 06241 xi_t[3] = getYS()->xi_t4( &EP_PredictorEPS ); 06242 06243 tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 06244 06245 hardMod_ = hardMod_ + hm.trace(); 06246 06247 Hh = Hh + (getPS()->d2Qodsdt4( &EP_PredictorEPS ))("ijmn") *(h_t[3])("mn") *Delta_lambda; 06248 06249 } 06250 06251 06252 06253 H = H + Hh; 06254 06255 06256 06257 06258 06259 // tensor temp3lower = dFods("ot")*R("otpq")*dQods("pq"); // Bug found, Z Cheng, Jan 2004 06260 06261 tensor temp3lower = dFods("ot")*R("otpq")*H("pq"); 06262 06263 temp3lower.null_indices(); 06264 06265 06266 06267 double lower = temp3lower.trace(); 06268 06269 06270 06271 lower = lower - hardMod_; // h 06272 06273 //opserr << " 2nd hardMod_ " << hardMod_ << "\n"; 06274 06275 06276 06277 06278 06279 //tensor upper = R("pqkl")*dQods("kl")*dFods("ij")*R("ijmn"); 06280 06281 // tensor upper11 = R("pqkl")*dQods("kl"); // Bug found, Z Cheng, Jan 2004 06282 06283 tensor upper11 = R("pqkl")*H("kl"); 06284 06285 tensor upper22 = dFods("ij")*R("ijmn"); 06286 06287 upper11.null_indices(); 06288 06289 upper22.null_indices(); 06290 06291 tensor upper = upper11("pq")*upper22("mn"); 06292 06293 06294 06295 upper.null_indices(); 06296 06297 tensor Ep = upper*(1./lower); 06298 06299 tensor Eep = R - Ep; // elastoplastic constitutive tensor 06300 06301 06302 06303 //opserr << "Eep :\n" << Eep; 06304 06305 06306 06307 //Set Elasto-Plastic stiffness tensor 06308 06309 EP_PredictorEPS.setEep(Eep); 06310 06311 EP_PredictorEPS.setConverged(TRUE); 06312 06313 06314 06315 //GZ out //set plastic strain and total strain 06316 06317 //GZ out incrPlasticStrain = dQods("kl")*Delta_lambda; 06318 06319 //GZ out incrPlasticStrain.null_indices(); 06320 06321 //GZ out 06322 06323 //GZ out PlasticStrain = startEPS.getStrain() - startEPS.getElasticStrain(); 06324 06325 //GZ out 06326 06327 //GZ out PlasticStrain = PlasticStrain + incrPlasticStrain; 06328 06329 //GZ out 06330 06331 //GZ out EP_PredictorEPS.setPlasticStrain(PlasticStrain); 06332 06333 //GZ out 06334 06335 //GZ out straintensor elastic_strain = strain_increment - incrPlasticStrain; // elastic strain increment 06336 06337 //GZ out straintensor estrain = startEPS.getElasticStrain(); //get old elastic strain 06338 06339 //GZ out //Guanzhou straintensor pstrain = EP_PredictorEPS.getPlasticStrain(); //get old plastic strain 06340 06341 //GZ out 06342 06343 //GZ out straintensor tstrain = startEPS.getStrain(); //get old total strain 06344 06345 //GZ out //Guanzhou pstrain = pstrain + PlasticStrain; 06346 06347 //GZ out estrain = estrain + elastic_strain; 06348 06349 //GZ out tstrain = tstrain + strain_increment; 06350 06351 //GZ out //opserr<< "Plastic: Total strain" << tstrain <<endln; 06352 06353 //GZ out 06354 06355 //GZ out //Setting de_p, de_e, total plastic, elastic strain, and total strain 06356 06357 //GZ out EP_PredictorEPS.setdPlasticStrain( incrPlasticStrain ); 06358 06359 //GZ out EP_PredictorEPS.setdElasticStrain( elastic_strain ); 06360 06361 //GZ out //EP_PredictorEPS.setPlasticStrain( pstrain ); 06362 06363 //GZ out EP_PredictorEPS.setElasticStrain( estrain ); 06364 06365 //GZ out EP_PredictorEPS.setStrain( tstrain ); 06366 06367 06368 06369 backwardEPS = EP_PredictorEPS; 06370 06371 06372 06373 //double f_backward = getYS()->f( &backwardEPS ); 06374 06375 //opserr << "\n************ Exit Backward = " << f_backward << "\n"; 06376 06377 06378 06379 } 06380 06381 06382 06383 //return Return_stress; 06384 06385 //Guanzhou out backwardEPS = EP_PredictorEPS; 06386 06387 06388 06389 //double p = (backwardEPS.getStress()).p_hydrostatic(); 06390 06391 //double ec = (backwardEPS.getec()) - (backwardEPS.getLam()) * log( p / (backwardEPS.getpo()) ); 06392 06393 //double st = (backwardEPS.getStrain()).Iinvariant1(); 06394 06395 //double pl_s = (backwardEPS.getPlasticStrain()).Iinvariant1(); 06396 06397 //double dpl_s = (backwardEPS.getdPlasticStrain()).Iinvariant1(); 06398 06399 //opserr << "P BE p=" << p << " ec " << ec << " e " << backwardEPS.gete() << " psi " << (backwardEPS.gete() - ec) << " strain " << st << " t_pl " << pl_s << " d_pl " << dpl_s << "\n"; 06400 06401 //opserr << "Backward EPState: \n" << backwardEPS; 06402 06403 return backwardEPS; 06404 06405 } 06406 06407 06408 06409 06410 06411 06412 06413 06414 06415 //================================================================================ 06416 06417 // New EPState using Forward Euler Subincrement Euler Algorithm 06418 06419 //================================================================================ 06420 06421 EPState Template3Dep::FESubIncrementation( const straintensor & strain_increment, 06422 06423 int number_of_subincrements) 06424 06425 { 06426 06427 06428 06429 EPState old_EPS( *(this->getEPS()) ); 06430 06431 EPState FESI_EPS( *(this->getEPS()) ); 06432 06433 //NDMaterial MP( material_point ); 06434 06435 //NDMaterial *MP = this->getCopy(); 06436 06437 //opserr << "in FESubIncrementation MP " << MP; 06438 06439 06440 06441 stresstensor back_stress; 06442 06443 stresstensor begin_stress = this->getEPS()->getStress(); 06444 06445 //stresstensor begin_stress = start_stress; 06446 06447 //::fprintf(stderr,"{"); 06448 06449 06450 06451 double sub = 1./( (double) number_of_subincrements ); 06452 06453 //stresstensor elastic_subincremental_stress = stress_increment * sub; 06454 06455 06456 06457 straintensor tempp = strain_increment; 06458 06459 straintensor elastic_subincremental_strain = tempp * sub; 06460 06461 straintensor total_strain = elastic_subincremental_strain; 06462 06463 //elastic_subincremental_stress.reportshort("SUB INCREMENT in stresses\n"); 06464 06465 //opserr << "INCREMENT strain " << strain_increment << endlnn ; 06466 06467 //opserr << "SUB INCREMENT strain " << elastic_subincremental_strain << endlnn ; 06468 06469 06470 06471 for( int steps=0 ; steps < number_of_subincrements ; steps++ ){ 06472 06473 06474 06475 //start_stress.reportshort("START stress\n"); 06476 06477 FESI_EPS = ForwardEulerEPState( elastic_subincremental_strain); 06478 06479 06480 06481 // Update the EPState in Template3Dep 06482 06483 this->setEPS( FESI_EPS ); 06484 06485 06486 06487 back_stress = FESI_EPS.getStress(); 06488 06489 //opserr.unsetf(ios::showpos); 06490 06491 //opserr << setw(4); 06492 06493 //opserr << "Step No. " << steps << " "; 06494 06495 06496 06497 //opserr.setPrecision(SCIENTIFIC); 06498 06499 //opserr.precision(3); 06500 06501 06502 06503 // opserr 06504 06505 // opserr.setf(ios::showpos); 06506 06507 // opserr.precision(3); 06508 06509 06510 06511 //opserr << setw(7); 06512 06513 //opserr << "p " << back_stress.p_hydrostatic() << " "; 06514 06515 //opserr << setw(7); 06516 06517 //opserr << "q " << back_stress.q_deviatoric() << " "; 06518 06519 //opserr << setw(7); 06520 06521 //opserr << " theta " << back_stress.theta() << " "; 06522 06523 //opserr << setw(7); 06524 06525 //opserr << "alfa1 " << FESI_EPS.getScalarVar(1) << " "; 06526 06527 //opserr << setw(7); 06528 06529 //opserr << "f = " << getYS()->f( &FESI_EPS ) << endlnn; 06530 06531 06532 06533 //begin_stress = back_stress; 06534 06535 //total_strain = total_strain + elastic_subincremental_strain; 06536 06537 06538 06539 } 06540 06541 06542 06543 // ::fprintf(stderr,"}"); 06544 06545 this->setEPS( old_EPS ); 06546 06547 return FESI_EPS; 06548 06549 06550 06551 } 06552 06553 06554 06555 06556 06557 06558 06559 //================================================================================ 06560 06561 // New EPState using Backward Euler Subincrement Euler Algorithm 06562 06563 //================================================================================ 06564 06565 EPState Template3Dep::BESubIncrementation( const straintensor & strain_increment, 06566 06567 int number_of_subincrements) 06568 06569 { 06570 06571 EPState old_EPS( *(this->getEPS()) ); 06572 06573 EPState BESI_EPS( *(this->getEPS()) ); 06574 06575 straintensor strain_incr = strain_increment; 06576 06577 //NDMaterial MP( material_point ); 06578 06579 //NDMaterial *MP = getCopy(); 06580 06581 //opserr << "in FESubIncrementation MP " << MP; 06582 06583 06584 06585 stresstensor back_stress; 06586 06587 stresstensor begin_stress = old_EPS.getStress(); 06588 06589 //stresstensor begin_stress = getEPS()->getStress(); 06590 06591 //stresstensor begin_stress = start_stress; 06592 06593 //::fprintf(stderr,"{"); 06594 06595 06596 06597 double sub = 1./( (double) number_of_subincrements ); 06598 06599 //stresstensor elastic_subincremental_stress = stress_increment * sub; 06600 06601 06602 06603 straintensor elastic_subincremental_strain = strain_incr * sub; 06604 06605 straintensor total_strain = elastic_subincremental_strain; 06606 06607 //elastic_subincremental_stress.reportshort("SUB INCREMENT in stresses\n"); 06608 06609 06610 06611 for( int steps=0 ; steps < number_of_subincrements ; steps++ ){ 06612 06613 06614 06615 //start_stress.reportshort("START stress\n"); 06616 06617 BESI_EPS = BackwardEulerEPState( elastic_subincremental_strain); 06618 06619 if ( BESI_EPS.getConverged() ) 06620 06621 this->setEPS( BESI_EPS ); 06622 06623 else { 06624 06625 //this->setEPS( BESI_EPS ); 06626 06627 opserr << "Template3Dep::BESubIncrementation failed to converge at " << steps << "th(of " 06628 06629 << number_of_subincrements << "step sub-BackwardEuler Algor.\n"; 06630 06631 //exit(1); 06632 06633 //g3ErrorHandler->fatal("Template3Dep::BESubIncrementation failed to converge using %d step sub-BackwardEuler Algor.", number_of_subincrements ); 06634 06635 //exit(1); 06636 06637 06638 06639 //back_stress = BESI_EPS.getStress(); 06640 06641 //opserr.unsetf(ios::showpos); 06642 06643 //opserr << setw(4); 06644 06645 //opserr << "Step No. " << steps << " "; 06646 06647 // 06648 06649 //opserr.setf(ios::scientific); 06650 06651 //opserr.setf(ios::showpos); 06652 06653 //opserr.precision(3); 06654 06655 //opserr << setw(7); 06656 06657 //opserr << " back-stress p " << back_stress.p_hydrostatic() << " "; 06658 06659 //opserr << setw(7); 06660 06661 //opserr << "q " << back_stress.q_deviatoric() << " "; 06662 06663 //opserr << setw(7); 06664 06665 //opserr << "alfa1 " << BESI_EPS.getScalarVar(1) << " "; 06666 06667 //opserr << setw(7); 06668 06669 //opserr << "f = " << MP->YS->f( &BESI_EPS ) << " "<< endlnn; 06670 06671 06672 06673 06674 06675 // opserr.setf(ios::scientific); 06676 06677 // opserr.setf(ios::showpos); 06678 06679 // opserr.precision(3); 06680 06681 opserr.width(7); 06682 06683 opserr << "\n intraStep: begin-stress p " << begin_stress.p_hydrostatic() << " "; 06684 06685 opserr.width(7); 06686 06687 opserr << "q " << begin_stress.q_deviatoric() << " "; 06688 06689 opserr.width(7); 06690 06691 opserr << "alfa1 " << old_EPS.getScalarVar(1) << " "; 06692 06693 //opserr.width(7); 06694 06695 //opserr << "f = " << MP->YS->f( &BESI_EPS ) << " "<< endlnn; 06696 06697 06698 06699 opserr << "strain increment " << strain_incr << endln; 06700 06701 06702 06703 //fflush(opserr); 06704 06705 break; 06706 06707 //exit(1); 06708 06709 //begin_stress = back_stress; 06710 06711 //total_strain = total_strain + elastic_subincremental_strain; 06712 06713 } 06714 06715 06716 06717 } 06718 06719 06720 06721 // ::fprintf(stderr,"}"); 06722 06723 //this->setEPS( old_EPS ); 06724 06725 return BESI_EPS; 06726 06727 06728 06729 } 06730 06731 06732 06733 06734 06735 06736 06738 06740 06742 06743 //tensor Template3Dep::ElasticStiffnessTensor( double E, double nu) const 06744 06745 // { 06746 06747 // tensor ret( 4, def_dim_4, 0.0 ); 06748 06749 // 06750 06751 // tensor I2("I", 2, def_dim_2); 06752 06753 // tensor I_ijkl = I2("ij")*I2("kl"); 06754 06755 // I_ijkl.null_indices(); 06756 06757 // tensor I_ikjl = I_ijkl.transpose0110(); 06758 06759 // tensor I_iljk = I_ijkl.transpose0111(); 06760 06761 // tensor I4s = (I_ikjl+I_iljk)*0.5; 06762 06763 // 06764 06765 // // Building elasticity tensor 06766 06767 // ret = (I_ijkl*((E*nu*2.0)/(2.0*(1.0+nu)*(1-2.0*nu)))) + (I4s*(E/((1.0+nu)))); 06768 06769 // 06770 06771 // return ret; 06772 06773 // } 06774 06775 // 06776 06777 // 06778 06780 06782 06784 06785 // 06786 06787 //tensor Template3Dep::ElasticComplianceTensor( double E, double nu) const 06788 06789 // { 06790 06791 // if (E == 0) { 06792 06793 // opserr << endln << "Ey = 0! Can't give you D!!" << endln; 06794 06795 // exit(1); 06796 06797 // } 06798 06799 // 06800 06801 // tensor ret( 4, def_dim_4, 0.0 ); 06802 06803 // //tensor ret; 06804 06805 // 06806 06807 // tensor I2("I", 2, def_dim_2); 06808 06809 // tensor I_ijkl = I2("ij")*I2("kl"); 06810 06811 // I_ijkl.null_indices(); 06812 06813 // tensor I_ikjl = I_ijkl.transpose0110(); 06814 06815 // tensor I_iljk = I_ijkl.transpose0111(); 06816 06817 // tensor I4s = (I_ikjl+I_iljk)*0.5; 06818 06819 // 06820 06821 // // Building compliance tensor 06822 06823 // ret = (I_ijkl * (-nu/E)) + (I4s * ((1.0+nu)/E)); 06824 06825 // 06826 06827 // return ret; 06828 06829 // } 06830 06831 // 06832 06833 06834 06835 //================================================================================ 06836 06837 // trying to find intersection point 06838 06839 // according to M. Crisfield's book 06840 06841 // "Non-linear Finite Element Analysis of Solids and Structures " 06842 06843 // chapter 6.6.1 page 168. 06844 06845 //================================================================================ 06846 06847 stresstensor Template3Dep::yield_surface_cross(const stresstensor & start_stress, 06848 06849 const stresstensor & end_stress) 06850 06851 { 06852 06853 // Bounds 06854 06855 double x1 = 0.0; 06856 06857 double x2 = 1.0; 06858 06859 06860 06861 // accuracy 06862 06863 double const TOL = 1.0e-9; 06864 06865 //opserr << "start_stress "<< start_stress; 06866 06867 //opserr << "end_stress " << end_stress; 06868 06869 //end_stress.reportshortpqtheta("end stress"); 06870 06871 06872 06873 double a = zbrentstress( start_stress, end_stress, x1, x2, TOL ); // Defined below 06874 06875 // ::printf("\n****\n a = %lf \n****\n",a); 06876 06877 06878 06879 stresstensor delta_stress = end_stress - start_stress; 06880 06881 stresstensor intersection_stress = start_stress + delta_stress * a; 06882 06883 //*** intersection_stress.reportshort("FINAL Intersection stress\n"); 06884 06885 06886 06887 return intersection_stress; 06888 06889 06890 06891 } 06892 06893 06894 06895 06896 06897 //================================================================================ 06898 06899 // Routine used by yield_surface_cross to 06900 06901 // find the stresstensor at cross point 06902 06903 //================================================================================ 06904 06905 06906 06907 double Template3Dep::zbrentstress(const stresstensor & start_stress, 06908 06909 const stresstensor & end_stress, 06910 06911 double x1, double x2, double tol) 06912 06913 { 06914 06915 double EPS = d_macheps(); 06916 06917 06918 06919 int iter; 06920 06921 double a=x1; 06922 06923 double b=x2; 06924 06925 double c=0.0; 06926 06927 double d=0.0; 06928 06929 double e=0.0; 06930 06931 double min1=0.0; 06932 06933 double min2=0.0; 06934 06935 double fa=func(start_stress, end_stress, a); 06936 06937 double fb=func(start_stress, end_stress, b); 06938 06939 //opserr << "fb = " << fb; 06940 06941 06942 06943 double fc=0.0; 06944 06945 double p=0.0; 06946 06947 double q=0.0; 06948 06949 double r=0.0; 06950 06951 double s=0.0; 06952 06953 double tol1=0.0; 06954 06955 double xm=0.0; 06956 06957 //::printf("\n############# zbrentstress iterations --\n"); 06958 06959 if (fb*fa > 0.0) 06960 06961 { 06962 06963 ::printf("\a\nRoot must be bracketed in ZBRENTstress\n"); 06964 06965 ::exit(1); 06966 06967 } 06968 06969 fc=fb; 06970 06971 for ( iter=1 ; iter<=ITMAX ; iter++ ) 06972 06973 { 06974 06975 //::printf("iter No. = %d ; b = %16.10lf\n", iter, b); 06976 06977 if (fb*fc > 0.0) 06978 06979 { 06980 06981 c=a; 06982 06983 fc=fa; 06984 06985 e=d=b-a; 06986 06987 } 06988 06989 if (fabs(fc) < fabs(fb)) 06990 06991 { 06992 06993 a=b; 06994 06995 b=c; 06996 06997 c=a; 06998 06999 fa=fb; 07000 07001 fb=fc; 07002 07003 fc=fa; 07004 07005 } 07006 07007 tol1=2.0*EPS*fabs(b)+0.5*tol; 07008 07009 xm=0.5*(c-b); 07010 07011 if (fabs(xm) <= tol1 || fb == 0.0) return b; 07012 07013 if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) 07014 07015 { 07016 07017 s=fb/fa; 07018 07019 if (a == c) 07020 07021 { 07022 07023 p=2.0*xm*s; 07024 07025 q=1.0-s; 07026 07027 } 07028 07029 else 07030 07031 { 07032 07033 q=fa/fc; 07034 07035 r=fb/fc; 07036 07037 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); 07038 07039 q=(q-1.0)*(r-1.0)*(s-1.0); 07040 07041 } 07042 07043 if (p > 0.0) q = -q; 07044 07045 p=fabs(p); 07046 07047 min1=3.0*xm*q-fabs(tol1*q); 07048 07049 min2=fabs(e*q); 07050 07051 if (2.0*p < (min1 < min2 ? min1 : min2)) 07052 07053 { 07054 07055 e=d; 07056 07057 d=p/q; 07058 07059 } 07060 07061 else 07062 07063 { 07064 07065 d=xm; 07066 07067 e=d; 07068 07069 } 07070 07071 } 07072 07073 else 07074 07075 { 07076 07077 d=xm; 07078 07079 e=d; 07080 07081 } 07082 07083 a=b; 07084 07085 fa=fb; 07086 07087 if (fabs(d) > tol1) 07088 07089 b += d; 07090 07091 else 07092 07093 b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1)); 07094 07095 fb=func(start_stress, end_stress, b); 07096 07097 } 07098 07099 //::printf("\a\nMaximum number of iterations exceeded in zbrentstress\n"); 07100 07101 return 0.0; // this just to full the compiler because of the warnings 07102 07103 } 07104 07105 07106 07107 07108 07109 //================================================================================ 07110 07111 // routine used by zbrentstress, takes an alfa and returns the 07112 07113 // yield function value for that alfa 07114 07115 //================================================================================ 07116 07117 double Template3Dep::func(const stresstensor & start_stress, 07118 07119 const stresstensor & end_stress, 07120 07121 double alfa ) 07122 07123 { 07124 07125 07126 07127 //EPState *tempEPS = getEPS()->newObj(); 07128 07129 EPState tempEPS( (*this->getEPS()) ); 07130 07131 stresstensor delta_stress = end_stress - start_stress; 07132 07133 stresstensor intersection_stress = start_stress + delta_stress*alfa; 07134 07135 07136 07137 tempEPS.setStress(intersection_stress); 07138 07139 07140 07141 //opserr << "*tempEPS" << *tempEPS; 07142 07143 07144 07145 double f = getYS()->f( &tempEPS ); 07146 07147 return f; 07148 07149 } 07150 07151 07152 07153 //================================================================================ 07154 07155 // Overloaded Insertion Operator 07156 07157 // prints an EPState's contents 07158 07159 //================================================================================ 07160 07161 OPS_Stream& operator<< (OPS_Stream& os, const Template3Dep & MP) 07162 07163 { 07164 07165 os << endln << "Template3Dep: " << endln; 07166 07167 os << "\ttag: " << MP.getTag() << endln; 07168 07169 os << "=================================================================" << endln; 07170 07171 MP.getYS()->print(); 07172 07173 MP.getPS()->print(); 07174 07175 MP.getEPS()->print(); 07176 07177 07178 07179 os << endln << "Scalar Evolution Laws: " << endln; 07180 07181 if ( MP.ELS1 ){ 07182 07183 os << "\nFor 1st scalar var:\n"; 07184 07185 MP.ELS1->print(); 07186 07187 } 07188 07189 07190 07191 if ( MP.ELS2 ){ 07192 07193 os << "\nFor 2nd scalar var:\n"; 07194 07195 MP.ELS2->print(); 07196 07197 } 07198 07199 07200 07201 if ( MP.ELS3 ){ 07202 07203 os << "\nFor 3rd scalar var:\n"; 07204 07205 MP.ELS3->print(); 07206 07207 } 07208 07209 07210 07211 if ( MP.ELS4 ){ 07212 07213 os << "\nFor 4th scalar var:\n"; 07214 07215 MP.ELS4->print(); 07216 07217 } 07218 07219 07220 07221 07222 07223 os << endln << "Tensorial Evolution Laws: " << endln; 07224 07225 if ( MP.ELT1 ){ 07226 07227 os << "\nFor 1st tensorial var:\n"; 07228 07229 MP.ELT1->print(); 07230 07231 } 07232 07233 if ( MP.ELT2 ){ 07234 07235 os << "\nFor 2nd tensorial var:\n"; 07236 07237 MP.ELT2->print(); 07238 07239 } 07240 07241 if ( MP.ELT3 ){ 07242 07243 os << "\nFor 3rd tensorial var:\n"; 07244 07245 MP.ELT3->print(); 07246 07247 } 07248 07249 if ( MP.ELT4 ){ 07250 07251 os << "\nFor 4th tensorial var:\n"; 07252 07253 MP.ELT4->print(); 07254 07255 } 07256 07257 07258 07259 os << endln; 07260 07261 return os; 07262 07263 } 07264 07265 07266 07267 07268 07269 //Guanzhou added, Apr2005 07270 07271 //EPState Template3Dep::NewBackwardEuler( const straintensor &strain_increment) 07272 07273 //{ 07274 07275 // EPState backwardEPS( *(this->getEPS()) ); 07276 07277 // EPState startEPS( *(this->getEPS()) ); 07278 07279 // 07280 07281 // double tol1 = 1e-8; 07282 07283 // double tol2 = 1e-8; 07284 07285 // 07286 07287 // int counter = 0; 07288 07289 // double Delta_lambda = 0.0; 07290 07291 // 07292 07293 // tensor E = ElasticStiffnessTensor(); 07294 07295 // 07296 07297 // straintensor strain_incr = strain_increment; 07298 07299 // straintensor total_strain = backwardEPS.getStrain(); 07300 07301 // total_strain = total_strain + strain_incr; 07302 07303 // backwardEPS.setStrain(total_strain);//watch out, iteration is carried out with constant total strain!!! 07304 07305 // 07306 07307 // straintensor plastic_strain = backwardEPS.getPlasticStrain(); 07308 07309 // straintensor trial_incr = total_strain - plastic_strain; 07310 07311 // 07312 07313 // stresstensor elastic_predictor = E("ijkl") * trial_incr("kl"); 07314 07315 // elastic_predictor.null_indices(); 07316 07317 // 07318 07319 // backwardEPS.setStress(elastic_predictor); 07320 07321 // 07322 07323 // tensor dFods(2, def_dim_2, 0.0); 07324 07325 // tensor dQods(2, def_dim_2, 0.0); 07326 07327 // tensor d2Qoverds2( 4, def_dim_4, 0.0); 07328 07329 // 07330 07331 // dQods = getPS()->dQods( &backwardEPS ); 07332 07333 // 07334 07335 // straintensor plastic_strain_residual = 07336 07337 // plastic_strain * (-1.0) + startEPS.getPlasticStrain() + dQods("kl") * Delta_lambda; 07338 07339 // 07340 07341 // //double hardMod_ = 0.0; 07342 07343 // Vector h_s(4); 07344 07345 // //double xi_s[4] = {0.0, 0.0, 0.0, 0.0}; 07346 07347 // stresstensor h_t[4]; 07348 07349 // //stresstensor xi_t[4]; 07350 07351 // 07352 07353 // //tensor Hh(2, def_dim_2, 0.0); 07354 07355 // //Of 1st scalar internal vars 07356 07357 // if ( getELS1() ) { 07358 07359 // h_s[0] = getELS1()->h_s( &backwardEPS, getPS()); 07360 07361 // //xi_s[0] = getYS()->xi_s1( &backwardEPS ); 07362 07363 // //hardMod_ = hardMod_ + h_s[0] * xi_s[0]; 07364 07365 // } 07366 07367 // 07368 07369 // //Of 2nd scalar internal vars 07370 07371 // if ( getELS2() ) { 07372 07373 // h_s[1] = getELS2()->h_s( &backwardEPS, getPS()); 07374 07375 // //xi_s[1] = getYS()->xi_s2( &backwardEPS ); 07376 07377 // //hardMod_ = hardMod_ + h_s[1] * xi_s[1]; 07378 07379 // } 07380 07381 // 07382 07383 // //Of 3rd scalar internal vars 07384 07385 // if ( getELS3() ) { 07386 07387 // h_s[2] = getELS3()->h_s( &backwardEPS, getPS()); 07388 07389 // //xi_s[2] = getYS()->xi_s3( &backwardEPS ); 07390 07391 // //hardMod_ = hardMod_ + h_s[2] * xi_s[2]; 07392 07393 // } 07394 07395 // 07396 07397 // //Of 4th scalar internal vars 07398 07399 // if ( getELS4() ) { 07400 07401 // h_s[3] = getELS4()->h_s( &backwardEPS, getPS()); 07402 07403 // //xi_s[3] = getYS()->xi_s4( &backwardEPS ); 07404 07405 // //hardMod_ = hardMod_ + h_s[3] * xi_s[3]; 07406 07407 // } 07408 07409 // 07410 07411 // //Of tensorial internal var 07412 07413 // // 1st tensorial var 07414 07415 // if ( getELT1() ) { 07416 07417 // h_t[0] = getELT1()->h_t( &backwardEPS, getPS()); 07418 07419 // //xi_t[0] = getYS()->xi_t1( &backwardEPS ); 07420 07421 // //tensor hm = (h_t[0])("ij") * (xi_t[0])("ij"); 07422 07423 // //hardMod_ = hardMod_ + hm.trace(); 07424 07425 // } 07426 07427 // 07428 07429 // // 2nd tensorial var 07430 07431 // if ( getELT2() ) { 07432 07433 // h_t[1] = getELT2()->h_t( &backwardEPS, getPS()); 07434 07435 // //xi_t[1] = getYS()->xi_t2( &backwardEPS ); 07436 07437 // //tensor hm = (h_t[1])("ij") * (xi_t[1])("ij"); 07438 07439 // //hardMod_ = hardMod_ + hm.trace(); 07440 07441 // } 07442 07443 // 07444 07445 // // 3rd tensorial var 07446 07447 // if ( getELT3() ) { 07448 07449 // h_t[2] = getELT3()->h_t( &backwardEPS, getPS()); 07450 07451 // //xi_t[2] = getYS()->xi_t3( &backwardEPS ); 07452 07453 // //tensor hm = (h_t[2])("ij") * (xi_t[2])("ij"); 07454 07455 // //hardMod_ = hardMod_ + hm.trace(); 07456 07457 // } 07458 07459 // 07460 07461 // // 4th tensorial var 07462 07463 // if ( getELT4() ) { 07464 07465 // h_t[3] = getELT4()->h_t( &backwardEPS, getPS()); 07466 07467 // //xi_t[3] = getYS()->xi_t4( &backwardEPS ); 07468 07469 // //tensor hm = (h_t[3])("ij") * (xi_t[3])("ij"); 07470 07471 // //hardMod_ = hardMod_ + hm.trace(); 07472 07473 // } 07474 07475 // 07476 07477 // int NS = backwardEPS.getNScalarVar(); 07478 07479 // int NT = backwardEPS.getNTensorVar(); 07480 07481 // 07482 07483 // Vector internal_variable_residual(NS+NT); 07484 07485 // 07486 07487 // stresstensor T; 07488 07489 // stresstensor new_T; 07490 07491 // stresstensor dT; 07492 07493 // 07494 07495 // tensor 07496 07497 // 07498 07499 // int ii; 07500 07501 // for (ii = 1; ii <= NS; ii++) { 07502 07503 // internal_variable_residual[ii-1] = 07504 07505 // backwardEPS.getScalarVar(ii)*(-1.0) + startEPS.getScalarVar(ii) + 07506 07507 // Delta_lambda * h_s[ii-1] ; 07508 07509 // } 07510 07511 // 07512 07513 // for (ii = 1; ii <= NT; ii++) { 07514 07515 // new_T = backwardEPS.getTensorVar(ii)*(-1.0); 07516 07517 // T = startEPS.getTensorVar(ii); 07518 07519 // dT = h_t[ii-1] * Delta_lambda; 07520 07521 // internal_variable_residual[NS+ii-1] = 07522 07523 // new_T.determinant() + T.determinant() + dT.determinant(); 07524 07525 // } 07526 07527 // 07528 07529 // double f = getYS()->f( &backwardEPS ); 07530 07531 // double residual = internal_variable_residual.Norm(); 07532 07533 // 07534 07535 // if ( (fabs(f) <= tol1) && residual <= tol2 ) { 07536 07537 // 07538 07539 // straintensor elastic_strain = backwardEPS.getElasticStrain(); 07540 07541 // elastic_strain = elastic_strain + strain_incr; 07542 07543 // backwardEPS.setElasticStrain(elastic_strain); 07544 07545 // backwardEPS.setEep(E); 07546 07547 // backwardEPS.setConverged(TRUE); 07548 07549 // //if ( getELT1() ) { 07550 07551 // // getELT1()->updateEeDm(&backwardEPS, -st_vol, 0.0); 07552 07553 // //} 07554 07555 // return backwardEPS; 07556 07557 // 07558 07559 // } else { 07560 07561 // 07562 07563 // while d2Qoverds2 = getPS()->d2Qods2( &backwardEPS ); 07564 07565 // 07566 07567 // return backwardEPS; 07568 07569 // } 07570 07571 // 07572 07573 //} 07574 07575 // 07576 07577 #endif 07578 07579 07580 07581 07582 07583 07584 07585 07586 07587 07588 07589 07590 07591 07592 07593 07594 07595 07596 07597 07598
Generated on Mon Oct 23 15:05:17 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents