YieldSurfaceSection2d.cppGo to the documentation of this file.00001 // @ rkaul@stanford.edu 00002 // @ ggd@stanford.edu 00003 00004 #include <math.h> 00005 00006 #include <YieldSurfaceSection2d.h> 00007 #include <Matrix.h> 00008 #include <Vector.h> 00009 #include <Channel.h> 00010 #include <FEM_ObjectBroker.h> 00011 #include <MatrixUtil.h> 00012 #include <stdlib.h> 00013 00014 #include <classTags.h> 00015 00016 ID YieldSurfaceSection2d::code(2); 00017 Vector YieldSurfaceSection2d::dele(2); 00018 Vector YieldSurfaceSection2d::surfaceForce(2); 00019 Matrix YieldSurfaceSection2d::G(2,1); 00020 Matrix YieldSurfaceSection2d::Ktp(2,2); 00021 00022 YieldSurfaceSection2d::YieldSurfaceSection2d(void) 00023 :SectionForceDeformation(0, SEC_TAG_YieldSurface2d), 00024 use_Kr_orig(true), ys(0), e(2), s(2), 00025 eCommit(2), sCommit(2), ks(2,2), 00026 use_Kr(true), split_step(false) 00027 { 00028 code(0) = SECTION_RESPONSE_P; // P is the first quantity 00029 code(1) = SECTION_RESPONSE_MZ; // Mz is the second 00030 } 00031 00032 YieldSurfaceSection2d::YieldSurfaceSection2d 00033 (int tag, int classtag, YieldSurface_BC *ptrys, bool use_kr) 00034 :SectionForceDeformation(tag, classtag), 00035 use_Kr_orig(use_kr), ys(0), e(2), s(2), 00036 eCommit(2), sCommit(2), ks(2,2), 00037 use_Kr(use_kr), split_step(false) 00038 { 00039 code(0) = SECTION_RESPONSE_P; // P is the first quantity 00040 code(1) = SECTION_RESPONSE_MZ; // Mz is the second 00041 00042 if(ptrys==0) 00043 { 00044 opserr << "WARNING - InelasticYS2DGNL(): ys1 = 0" << endln; 00045 } 00046 else 00047 { 00048 ys = ptrys->getCopy(); 00049 ys->setTransformation(1, 0, 1, -1); // x-axis is Mz, y-axis is P 00050 // ys->setEleInfo(getTag(), 1); 00051 } 00052 } 00053 00054 00055 YieldSurfaceSection2d::~YieldSurfaceSection2d(void) 00056 { 00057 if (ys != 0) 00058 delete ys; 00059 } 00060 00061 int 00062 YieldSurfaceSection2d::commitState(void) 00063 { 00064 eCommit = e; 00065 sCommit = s; 00066 ys->commitState(s); 00067 split_step = false; 00068 00069 return 0; 00070 } 00071 00072 int 00073 YieldSurfaceSection2d::revertToLastCommit(void) 00074 { 00075 e = eCommit; 00076 s = sCommit; 00077 ys->revertToLastCommit(); 00078 00079 return 0; 00080 } 00081 00082 int 00083 YieldSurfaceSection2d::revertToStart(void) 00084 { 00085 eCommit.Zero(); 00086 sCommit.Zero(); 00087 // ys->revertToStart(); 00088 00089 return 0; 00090 } 00091 00092 int 00093 YieldSurfaceSection2d::setTrialSectionDeformation (const Vector &def) 00094 { 00095 ys->update(); // important 00096 use_Kr = use_Kr_orig; 00097 // split_step = false; 00098 // once determined, leave it till convergence 00099 00100 e = def; 00101 dele = e - eCommit; 00102 00103 getSectionStiffness(ks); 00104 double EA = ks(0,0); 00105 double EI = ks(1,1); 00106 00107 s(0) = sCommit(0) + EA*dele(0); 00108 s(1) = sCommit(1) + EI*dele(1); 00109 00110 if(ys->getTrialForceLocation(s) <= 0) 00111 return 0; 00112 00113 00114 // case: if it shoots through 00115 // use dF return to surface 00116 int driftOld = ys->getCommitForceLocation(); 00117 00118 if(driftOld < 0) // Inside 00119 { 00120 use_Kr = false; 00121 split_step = true; 00122 00123 surfaceForce = s; 00124 ys->setToSurface(surfaceForce, ys->dFReturn); //dFReturn, ConstantYReturn, RadialReturn 00125 ys->getTrialGradient(G, surfaceForce); 00126 } 00127 // Now we know that force point has drifted from the surface 00128 else if(driftOld == 0) // On surface 00129 { 00130 ys->getCommitGradient(G); 00131 surfaceForce = sCommit; 00132 } 00133 else // driftOld outside - bug or bad constraints or continued from not converged state 00134 { 00135 opserr << "WARNING: YieldSurfaceSection2d::setTrialSectionDeformation, driftOld outside [" << getTag()<<"]\n"; 00136 } 00137 00138 double dF_in0 = s(0) - surfaceForce(0); 00139 double dF_in1 = s(1) - surfaceForce(1); 00140 00141 double g0 = G(0,0); 00142 double g1 = G(1,0); 00143 00144 Ktp(0,0) = EA; 00145 Ktp(1,1) = EI; 00146 ys->addPlasticStiffness(Ktp); 00147 00148 double inv = 1/(Ktp(0,0)*g0*g0 + Ktp(1,1)*g1*g1); 00149 00150 double lamda = inv*(g0*dF_in0 + g1*dF_in1); 00151 if(fabs(lamda) < 1e-8) lamda = 0.0; // to get rid of -1e-15 etc 00152 00153 if(lamda < 0) 00154 { 00155 use_Kr = false; 00156 lamda = 0.0; 00157 } 00158 00159 int grow = ys->modifySurface(lamda, surfaceForce, G); 00160 00161 // used to do: (not tested shrinking yet) 00162 // if(grow < 0) 00163 // algo = ys->ConstantYReturn; 00164 // else 00165 // algo = algo_orig; 00166 00167 if(use_Kr) 00168 { 00169 ks(0,0) = EA - EA*EA*g0*g0*inv; 00170 ks(0,1) = -1*EA*g0*g1*EI*inv; 00171 ks(1,0) = ks(0,1); 00172 ks(1,1) = EI - EI*EI*g1*g1*inv; 00173 } 00174 if(split_step) 00175 { 00176 s(0) = s(0) - EA*g0*lamda; 00177 s(1) = s(1) - EI*g1*lamda; 00178 } 00179 else 00180 { 00181 if(use_Kr) 00182 s = surfaceForce + ks*dele; 00183 } 00184 00185 ys->setToSurface(s, ys->ConstantYReturn); 00186 // used to do centroid return 00187 // then force-balance using ConstantYReturn 00188 // comp/tension issue: always use constantP 00189 00190 return 0; 00191 } 00192 00193 const Vector & 00194 YieldSurfaceSection2d::getSectionDeformation (void) 00195 { 00196 return e; 00197 } 00198 00199 const Vector & 00200 YieldSurfaceSection2d::getStressResultant (void) 00201 { 00202 return s; 00203 } 00204 00205 const Matrix & 00206 YieldSurfaceSection2d::getSectionTangent(void) 00207 { 00208 return ks; 00209 } 00210 00211 const Matrix & 00212 YieldSurfaceSection2d::getSectionFlexibility (void) 00213 { 00214 static Matrix fs(2,2); 00215 invert2by2Matrix(ks, fs); 00216 00217 return fs; 00218 } 00219 00220 const ID& 00221 YieldSurfaceSection2d::getType () 00222 { 00223 return code; 00224 } 00225 00226 int 00227 YieldSurfaceSection2d::getOrder () const 00228 { 00229 return 2; 00230 } 00231 00232 int 00233 YieldSurfaceSection2d::sendSelf(int commitTag, Channel &theChannel) 00234 { 00235 return -1; 00236 } 00237 00238 int 00239 YieldSurfaceSection2d::recvSelf(int commitTag, Channel &theChannel, 00240 FEM_ObjectBroker &theBroker) 00241 { 00242 return -1; 00243 } 00244 00245 void 00246 YieldSurfaceSection2d::Print(OPS_Stream &s, int flag) 00247 { 00248 s << "YieldSurfaceSection2d, tag: " << this->getTag() << endln; 00249 s << "\tYield Surface:" << *ys << endln; 00250 s << "\tSection Force:" << sCommit; 00251 s << "\tSection Defom:" << eCommit; 00252 }
Generated on Mon Oct 23 15:05:19 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents