Cosseratstresst.cppGo to the documentation of this file.00001 00002 //################################################################################ 00003 //# COPY-YES (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: stress tensor with all necessery functions # 00006 //# CLASS: stresstensor # 00007 //# # 00008 //# VERSION: # 00009 //# LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) # 00010 //# TARGET OS: DOS || UNIX || . . . # 00011 //# DESIGNER(S): Alireza Tabarrei # 00012 //# PROGRAMMER(S): Alireza Tabarrei # 00013 //# # 00014 //# # 00015 //# DATE: June 2004 # 00016 //# UPDATE HISTORY: # 00017 //# # 00018 //# # 00019 //# # 00020 //# # 00021 //################################################################################ 00022 //*/ 00023 00024 00025 #ifndef COSSERATSTRESSTENSOR_CPP 00026 #define COSSERATSTRESSTENSOR_CPP 00027 00028 #include "Cosseratstresst.h" 00029 00030 #include <iomanip> 00031 using std::ios; 00032 00033 // just send appropriate arguments to the base constructor 00034 //############################################################################## 00035 Cosseratstresstensor::Cosseratstresstensor (int rank_of_tensor, double initval): 00036 BJtensor(rank_of_tensor, Cosserat_def_dim_2, initval) { } // default constructor 00037 00038 00039 //############################################################################## 00040 Cosseratstresstensor::Cosseratstresstensor ( double *values ): 00041 BJtensor( 2, Cosserat_def_dim_2, values) { } 00042 00043 //############################################################################## 00044 Cosseratstresstensor::Cosseratstresstensor ( double initvalue ): 00045 BJtensor( 2, Cosserat_def_dim_2, initvalue) { } 00046 00047 //############################################################################## 00048 Cosseratstresstensor::Cosseratstresstensor( const Cosseratstresstensor & x ): 00049 BJtensor("NO") 00050 { 00051 x.pc_nDarray_rep->n++; // tell the rval it has another reference 00052 // x.reference_count(+1); // we're adding another reference. 00053 pc_nDarray_rep = x.pc_nDarray_rep; // point to the new tensor_rep. 00054 // add the indices 00055 indices1 = x.indices1; 00056 indices2 = x.indices2; 00057 } 00058 00059 00060 //############################################################################## 00061 Cosseratstresstensor::Cosseratstresstensor( const tensor & x): 00062 BJtensor( x ) 00063 { } // copy-initializer 00064 00065 //############################################################################## 00066 Cosseratstresstensor::Cosseratstresstensor( const nDarray & x): 00067 BJtensor ( x ) 00068 { } // copy-initializer 00069 00070 00071 00072 // IT IS NOT INHERITED so must be defined in all derived classes 00073 // See ARM page 277. 00074 //############################################################################## 00075 // Cosseratstresstensor::~Cosseratstresstensor() 00076 // { 00077 // if (reference_count(-1) <= 0) // Zhaohui changed == 0 to <= 0 // if reference count goes to 0 00078 // { 00079 // // DEallocate memory of the actual nDarray 00080 // // delete [pc_nDarray_rep->pc_nDarray_rep->total_numb] pc_nDarray_rep->pd_nDdata; 00081 // // nema potrebe za brojem clanova koji se brisu## see ELLIS & STROUSTRUP $18.3 00082 // // and note on the p.65($5.3.4) 00083 // delete [] data(); 00084 // delete [] dim(); 00085 // delete pc_nDarray_rep; 00086 // } 00087 // } 00088 // 00089 00090 //############################################################################## 00091 // IT IS NOT INHERITED so must be defined in all derived classes 00092 // See ARM page 306. 00093 //############################################################################## 00094 Cosseratstresstensor Cosseratstresstensor::operator=( const Cosseratstresstensor & rval) 00095 { 00096 rval.pc_nDarray_rep->n++; // tell the rval it has another reference 00097 // rval.reference_count(+1); // tell the rval it has another reference 00098 // /* It is important to increment the reference_counter in the new 00099 // tensor before decrementing the reference_counter in the 00100 // old tensor_rep to ensure proper operation when assigning a 00101 // tensor_rep to itself ( after ARKoenig JOOP May/June '90 ) */ 00102 // clean up current value; 00103 if( reference_count(-1) == 0) // if nobody else is referencing us. 00104 { 00105 // DEallocate memory of the actual tensor 00106 delete [] data(); 00107 delete [] dim(); 00108 delete pc_nDarray_rep; 00109 } 00110 // connect to new value 00111 pc_nDarray_rep = rval.pc_nDarray_rep; // point at the rval tensor_rep 00112 // Temporary out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00113 // null indices in the rval 00114 // rval.indices1 = NULL; 00115 // rval.indices2 = NULL; 00116 // rval.null_indices(); 00117 // Temporary out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00118 this->null_indices(); 00119 return *this; 00120 } 00121 00122 // IT IS NOT INHERITED so must be defined in all derived classes 00123 // See ARM page 306. 00124 //############################################################################## 00125 Cosseratstresstensor Cosseratstresstensor::operator=( const tensor & rval) 00126 { 00127 rval.pc_nDarray_rep->n++; // tell the rval it has another reference 00128 // rval.pc_nDarray_rep->n++; // tell the rval it has another reference 00129 // rval.reference_count(+1); // tell the rval it has another reference 00130 // /* It is important to increment the reference_counter in the new 00131 // tensor before decrementing the reference_counter in the 00132 // old tensor_rep to ensure proper operation when assigning a 00133 // tensor_rep to itself ( after ARKoenig JOOP May/June '90 ) */ 00134 00135 // clean up current value; 00136 // if(--pc_nDarray_rep->n == 0) // if nobody else is referencing us. 00137 if( reference_count(-1) == 0) // if nobody else is referencing us. 00138 { 00139 // DEallocate memory of the actual tensor 00140 // delete [pc_tensor_rep->pc_tensor_rep->total_numb] pc_tensor_rep->pd_nDdata; 00141 // nema potrebe za brojem clanova koji se brisu## see ELLIS & STROUSTRUP $18.3 00142 // and note on the p.65($5.3.4) 00143 // delete pc_nDarray_rep->pd_nDdata; 00144 delete [] data(); 00145 delete [] dim(); 00146 // ovo ne smem da brisem jer nije dinamicki alocirano 00147 // delete pc_tensor_rep->indices; 00148 delete pc_nDarray_rep; 00149 } 00150 00151 // connect to new value 00152 pc_nDarray_rep = rval.pc_nDarray_rep; // point at the rval tensor_rep 00153 // Temporary out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00154 // null indices in the rval 00155 // rval.indices1 = NULL; 00156 // rval.indices2 = NULL; 00157 // rval.null_indices(); 00158 // Temporary out !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00159 this->null_indices(); 00160 return *this; 00161 } 00162 00163 00164 // IT IS NOT INHERITED so must be defined in all derived classes 00165 // See ARM page 306. 00166 //############################################################################## 00167 Cosseratstresstensor Cosseratstresstensor::operator=( const nDarray & rval) 00168 { 00169 rval.pc_nDarray_rep->n++; // tell the rval it has another reference 00170 // rval.reference_count(+1); // tell the rval it has another reference 00171 // /* It is important to increment the reference_counter in the new 00172 // tensor before decrementing the reference_counter in the 00173 // old tensor_rep to ensure proper operation when assigning a 00174 // tensor_rep to itself ( after ARKoenig JOOP May/June '90 ) */ 00175 00176 if( reference_count(-1) == 0) // if nobody else is referencing us. 00177 { 00178 delete [] data(); 00179 delete [] dim(); 00180 delete pc_nDarray_rep; 00181 } 00182 00183 // connect to new value 00184 pc_nDarray_rep = rval.pc_nDarray_rep; // point at the rval tensor_rep 00185 return *this; 00186 } 00187 00188 00189 //############################################################################## 00190 // makes a complete new copy of Cosseratstresstensor!! 00191 Cosseratstresstensor Cosseratstresstensor::deep_copy(void) 00192 { 00193 return Cosseratstresstensor(this->data()); // call constructor and return it ! 00194 } 00195 //..//############################################################################## 00196 //../ returns a pointer to the Cosseratstresstensor!! 00197 //..Cosseratstresstensor * Cosseratstresstensor::p_deep_copy(void) 00198 //.. { 00199 //.. Cosseratstresstensor temp = this->deep_copy(); 00200 //.. return &temp; // call constructor and return it ! 00201 //.. } 00202 00203 //ini //############################################################################## 00204 //ini // use "from" and initialize already allocated stress tensor from "from" values 00205 //ini void Cosseratstresstensor::Initialize( const Cosseratstresstensor & from ) 00206 //ini { 00207 //ini // copy only data because everything else is default 00208 //ini for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00209 //ini this->pc_nDarray_rep->pd_nDdata[i] = from.pc_nDarray_rep->pd_nDdata[i] ; 00210 //ini } 00211 00212 //___//############################################################################## 00213 //___//############################################################################## 00214 //___//############################################################################## 00215 //___Cosseratstresstensor & Cosseratstresstensor::operator()(short ir, short is, short it, 00216 //___ short tr, short ts, short tt ) 00217 //___// another overloading of operator() . . . // overloaded for THREE arguments 00218 //___ { 00219 //___ short where = ir - 1; 00220 //___ where = where*ts + is - 1; 00221 //___ where = where*tt + it - 1; 00222 //___ 00223 //___//::printf(" w=%ld ",where); 00224 //___ Cosseratstresstensor *p_value = this + where; 00225 //___ return (*p_value); 00226 //___ } 00227 00228 //############################################################################## 00229 // invariants of the stress tensor // Chen W.F. "plasticity for 00230 double Cosseratstresstensor::Iinvariant1() const // Structural Engineers" 00231 { 00232 return (cval(1,1)+cval(2,2)+cval(3,3) + cval(4,4)+cval(5,5)+cval(6,6) ); 00233 } 00234 00235 00237 //############################################################################## 00238 double Cosseratstresstensor::Iinvariant2() const 00239 { 00240 return (cval(2,2)*cval(3,3)-cval(3,2)*cval(2,3)+ 00241 cval(1,1)*cval(3,3)-cval(3,1)*cval(1,3)+ 00242 cval(1,1)*cval(2,2)-cval(2,1)*cval(1,2)); 00243 } 00244 00245 //############################################################################## 00246 double Cosseratstresstensor::Iinvariant3() const 00247 { 00248 00249 double I3 = cval(1,1)*cval(2,2)*cval(3,3) + 00250 cval(1,2)*cval(2,3)*cval(3,1) + 00251 cval(1,3)*cval(2,1)*cval(3,2) - 00252 cval(1,3)*cval(2,2)*cval(3,1) - 00253 cval(1,2)*cval(2,1)*cval(3,3) - 00254 cval(1,1)*cval(2,3)*cval(3,2) ; 00255 00256 return I3; 00257 // return ( this->determinant()); 00258 } 00259 00260 00261 //############################################################################## 00262 // invariants of the deviatoric stress tensor 00263 double Cosseratstresstensor::Jinvariant1() const 00264 { 00265 return (0.0); 00266 } 00267 00268 //############################################################################## 00269 double Cosseratstresstensor::Jinvariant2() const 00270 { 00271 double temp1 = (Iinvariant1()*Iinvariant1()-3.0*Iinvariant2())*ONEOVERTHREE; 00272 // double EPS = sqrt(d_macheps()); 00273 // if ( temp1 < 0.0 || fabs(temp1) < EPS ) 00274 if ( temp1 < 0.0 ) 00275 { // this is because it might be close 00276 temp1 = 0.0; // to zero ( -1e-19 ) numericaly 00277 } // but sqrt() does not accept it 00278 return ( temp1 ); // as (-) and theoreticaly J2d>0 00279 } 00280 00281 //############################################################################## 00282 double Cosseratstresstensor::Jinvariant3() const 00283 { 00284 double I1 = Iinvariant1(); 00285 double I2 = Iinvariant2(); 00286 double I3 = Iinvariant3(); 00287 00288 return ( (2.0*I1*I1*I1 - 9.0*I1*I2 + 27.0*I3)/27.0 ); 00289 } 00290 00291 //############################################################################## 00292 Cosseratstresstensor Cosseratstresstensor::principal() const 00293 { 00294 Cosseratstresstensor ret; 00295 00296 double p_ = this->p_hydrostatic(); 00297 double q_ = this->q_deviatoric(); 00298 double theta_ = this->theta(); 00299 //old while ( theta_ >= 2.0*PI ) 00300 //old theta_ = theta_ - 2.0*PI; // if bigger than full cycle 00301 //old while ( theta_ >= 4.0*PI/3.0 ) 00302 //old theta_ = theta_ - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00303 //old while ( theta_ >= 2.0*PI/3.0 ) 00304 //old theta_ = theta_ - 2.0*PI/3.0; // if bigger than two third of half cycle 00305 //old while ( theta_ >= PI/3.0 ) 00306 //old theta_ = 2.0*PI/3.0 - theta_; // if bigger than one third of half cycle 00307 00308 // ret.report("ret"); 00309 // double sqrtJ2D = q/3.0; 00310 // double 2osqrt3 = 2.0/sqrt(3.0); 00311 double temp = q_*TWOOVERTHREE; 00312 00313 double ct = cos( theta_ ); 00314 double ctm = cos( theta_ - TWOOVERTHREE*PI ); 00315 double ctp = cos( theta_ + TWOOVERTHREE*PI ); 00316 00317 ret.val(1,1) = - p_ + temp*ct; // - becasue p is p = -1/3 sigma_ij delta_ii 00318 ret.val(2,2) = - p_ + temp*ctm; 00319 ret.val(3,3) = - p_ + temp*ctp; 00320 00321 // ret.report("ret"); 00322 return ret; 00323 00324 //.. static complex ac[4]; // Chen W.F. "plasticity for 00325 //.. static complex roots[4]; // Structural Engineers" 00326 //.. int polish = 1 ; // page 53 00327 //.. int m = 3; 00328 //.. 00329 //.. ac[0] = complex( -(this->Iinvariant3()), 0.0 ); 00330 //.. ac[1] = complex( (this->Iinvariant2()), 0.0); 00331 //.. ac[2] = complex( -(this->Iinvariant1()), 0.0); 00332 //.. ac[3] = complex( 1.0, 0.0); 00333 //.. 00334 //..// what was obtained for coefficients 00335 //..//DEBUGprint ::printf("ac[0].r = %lf ac[0].i = %lf\n", real(ac[0]), imag(ac[0])); 00336 //..//DEBUGprint ::printf("ac[1].r = %lf ac[1].i = %lf\n", real(ac[1]), imag(ac[1])); 00337 //..//DEBUGprint ::printf("ac[2].r = %lf ac[2].i = %lf\n", real(ac[2]), imag(ac[2])); 00338 //..//DEBUGprint ::printf("ac[3].r = %lf ac[3].i = %lf\n", real(ac[3]), imag(ac[3])); 00339 //..//DEBUGprint ::printf("m = %d\n",m); 00340 //.. 00341 //.. zroots( ac, m, roots, polish); 00342 //.. 00343 //..//DEBUGprint ::printf("\nroots[1].r=%lf roots[1].i=%lf\n",real(roots[1]),imag(roots[1])); 00344 //..//DEBUGprint ::printf("roots[2].r=%lf roots[2].i=%lf\n",real(roots[2]),imag(roots[2])); 00345 //..//DEBUGprint ::printf("roots[3].r=%lf roots[3].i=%lf\n",real(roots[3]),imag(roots[3])); 00346 //.. 00347 //.. 00348 //.. Cosseratstresstensor principal(0.0); 00349 //.. 00350 //.. principal.val(1,1) = real(roots[3]); // since they are sorted by 00351 //.. principal.val(2,2) = real(roots[2]); // the zroot function in ascending 00352 //.. principal.val(3,3) = real(roots[1]); // order . . . 00353 //.. // sig1>sig2>sig3 00354 //.. return principal; 00355 //.. 00356 //.. 00357 } 00358 00359 //############################################################################## 00360 Cosseratstresstensor Cosseratstresstensor::deviator() const 00361 { 00362 BJtensor I2("I", 2, def_dim_2); // Kronecker delta \delta_{ij} 00363 Cosseratstresstensor st_vol = I2 * (this->trace()*(0.333333333)); 00364 Cosseratstresstensor st_dev = (*this) - st_vol; 00365 00366 return st_dev; 00367 } 00368 00369 00370 00371 //############################################################################## 00372 double Cosseratstresstensor::sigma_octahedral() const // Chen W.F. "plasticity for 00373 { // Structural Engineers" 00374 return ( this->Iinvariant1()*ONEOVERTHREE ); // page 59-60 00375 } 00376 00377 //############################################################################## 00378 double Cosseratstresstensor::tau_octahedral() const // Chen W.F. "plasticity for 00379 { // Structural Engineers" 00380 return(sqrt(TWOOVERTHREE*(this->Jinvariant2())));// page 59-60 00381 } 00382 00383 00384 //############################################################################## 00385 double Cosseratstresstensor::ksi() const // Chen W.F. "plasticity for 00386 { // Structural Engineers" 00387 return( (this->Iinvariant1())/sqrt(3.0) ); // page 66 00388 } 00389 00390 //############################################################################## 00391 double Cosseratstresstensor::xi() const // Chen W.F. "plasticity for 00392 { // Structural Engineers" 00393 return( (this->Iinvariant1())/sqrt(3.0) ); // page 66 00394 } 00395 00396 //############################################################################## 00397 double Cosseratstresstensor::ro() const // Chen W.F. "plasticity for 00398 { // Structural Engineers" 00399 double temp1 = this->Jinvariant2(); // page 68 00400 double EPS = pow(d_macheps(),(1./2.)); 00401 if ( temp1 < 0.0 && fabs(temp1) < EPS ) 00402 { 00403 temp1 = 0.0; 00404 } 00405 return( sqrt(2.0*(temp1))); 00406 } 00407 //############################################################################## 00408 double Cosseratstresstensor::rho() const // Chen W.F. "plasticity for 00409 { // Structural Engineers" 00410 double temp1 = this->Jinvariant2(); // page 68 00411 double EPS = pow(d_macheps(),(1./2.)); 00412 if ( temp1 < 0.0 && fabs(temp1) < EPS ) 00413 { 00414 temp1 = 0.0; 00415 } 00416 return( sqrt(2.0*(temp1))); 00417 } 00418 00419 //############################################################################## 00420 double Cosseratstresstensor::p_hydrostatic() const // Desai "Constitutive Laws 00421 { // for Engineering Materials" 00422 return( - (this->Iinvariant1())*ONEOVERTHREE ); // page 283 00423 } //sign (-) because tension is positive 00424 00425 00426 //############################################################################## 00427 double Cosseratstresstensor::q_deviatoric() const // Desai "Constitutive Laws 00428 { // for Engineering Materials" 00429 double temp1 = this->Jinvariant2(); // page 283 00430 // double EPS = pow(d_macheps(),(1./2.)); 00431 // if ( temp1 < 0.0 || fabs(temp1) < EPS ) 00432 if ( temp1 < 0.0 ) 00433 { 00434 temp1 = 0.0; 00435 } 00436 00437 return( sqrt(3.0*temp1) ); 00438 } 00439 00440 //############################################################################## 00441 double Cosseratstresstensor::theta( ) const // Chen W.F. "plasticity for 00442 { // Structural Engineers" 00443 // page 70 00444 double EPS = pow(d_macheps(),(1./4.)); 00445 double temp1 = (3.0*sqrt(3.0)/2.0); 00446 double J3D = (this->Jinvariant3()); 00447 double J2D = (this->Jinvariant2()); 00448 if ( J2D < 0.0 || fabs(J2D) < EPS ) 00449 { 00450 J2D = 0.0; 00451 } 00452 double temp4 = J2D * J2D * J2D; 00453 double temp5 = temp1 * J3D; 00454 double temp6 = sqrt(temp4); 00455 if ( (fabs(temp6)) <= fabs(temp5 * EPS) ) return ( 0.000001 );// slight perturbation because of 1/0 KERU06sep96 00456 if ( (fabs(temp6)) < EPS ) return ( 0.000001 );// slight perturbation because of 1/0 KERU06sep96 00457 double temp7 = temp5 / temp6; 00458 double tempabs1 = (fabs(temp7-1.0)); 00459 double tempabs2 = (fabs(temp7+1.0)); 00460 if ( tempabs1 < 0.1 ) return ( 0.000001/3.0 );// slight perturbation because of 1/0 KERU06sep96 00461 // if ( tempabs2 < 0.001 ) return ( PI/3.0 ); 00462 if ( tempabs2 < 0.1 ) return ( 3.14159/3.0 );// slight perturbation because of 1/0 KERU06sep96 00463 if ( temp7>1.0 || temp7<-1.0 ) 00464 { 00465 ::fprintf(stderr,"\n something is wrong in Cosseratstresstensor::theta() (temp7>1.0||temp7<-1.0)\n"); 00466 ::fprintf(stderr,"returning Pi/3.\n"); 00467 00468 ::fprintf(stderr,"temp1 = %.20e\n", temp1); 00469 ::fprintf(stderr,"J3D = %.20e\n", J3D); 00470 ::fprintf(stderr,"J2D = %.20e\n", J2D); 00471 ::fprintf(stderr,"temp4 = %.20e\n", temp4); 00472 ::fprintf(stderr,"tempabs1 = %.20e\n", tempabs1); 00473 ::fprintf(stderr,"tempabs2 = %.20e\n", tempabs2); 00474 ::fprintf(stderr,"temp5 = %.20e\n", temp5); 00475 ::fprintf(stderr,"temp6 = %.20e\n", temp6); 00476 ::fprintf(stderr,"temp7 = %.20e\n", temp7); 00477 00478 ::fprintf(stdout,"\n something is wrong in Cosseratstresstensor::theta() (temp7>1.0||temp7<-1.0)\n"); 00479 ::fprintf(stdout,"returning Pi/3.\n"); 00480 this->print("s","Cosseratstresstensor s"); 00481 ::fprintf(stdout,"temp1 = %.20e\n", temp1); 00482 ::fprintf(stdout,"J3D = %.20e\n", J3D); 00483 ::fprintf(stdout,"J2D = %.20e\n", J2D); 00484 ::fprintf(stdout,"temp4 = %.20e\n", temp4); 00485 ::fprintf(stdout,"tempabs1 = %.20e\n", tempabs1); 00486 ::fprintf(stdout,"tempabs2 = %.20e\n", tempabs2); 00487 ::fprintf(stdout,"temp5 = %.20e\n", temp5); 00488 ::fprintf(stdout,"temp6 = %.20e\n", temp6); 00489 ::fprintf(stdout,"temp7 = %.20e\n", temp7); 00490 00491 return ( 3.14159/ 3.0 );// slight perturbation because of 1/0 KERU06sep96 00492 // exit(1); 00493 } 00494 double temp8 = acos(temp7); 00495 double temp9 = temp8 * ONEOVERTHREE; 00496 00497 return ( temp9 ); 00498 } 00499 00500 //############################################################################## 00501 double Cosseratstresstensor::thetaPI() const 00502 { 00503 double thetaPI = theta() / PI; 00504 return thetaPI; 00505 } 00506 00507 00511 // double Cosseratstresstensor::func( Cosseratstresstensor & start_stress, 00512 // Cosseratstresstensor & end_stress, 00513 // Material_Model & YC, 00514 // double alfa ) 00515 // { 00516 // Cosseratstresstensor delta_stress = end_stress - start_stress; 00517 // Cosseratstresstensor intersection_stress = start_stress + delta_stress*alfa; 00518 // double f = YC.f(intersection_stress); 00519 // return f; 00520 // } 00521 00522 00523 00524 00525 //############################################################################# 00526 tensor Cosseratstresstensor::dpoverds( void ) const 00527 { 00528 tensor ret(2, def_dim_2, 0.0); 00529 tensor I2("I", 2, def_dim_2); 00530 ret = I2*(-1.0/3.0); 00531 ret.null_indices(); 00532 return ret; 00533 } 00534 00535 //############################################################################# 00536 tensor Cosseratstresstensor::dqoverds( void ) const 00537 { 00538 00539 // Cosseratstresstensor stress = EPS->getStress(); 00540 00541 tensor ret(2, def_dim_2, 0.0); 00542 double q = this->q_deviatoric(); 00543 Cosseratstresstensor s( 0.0); 00544 //... double J2D = stress.Jinvariant2(); 00545 //... double temp1 = sqrt(J2D); 00546 //... double temp2 = sqrt(3.0)/(2.0*temp1); 00547 double temp2 = (3.0/2.0)*(1/q); 00548 s = this->deviator(); 00549 ret = s*temp2; 00550 ret.null_indices(); 00551 return ret; 00552 } 00553 00554 //############################################################################# 00555 tensor Cosseratstresstensor::dthetaoverds( void ) const 00556 { 00557 // Cosseratstresstensor stress = EPS->getStress(); 00558 00559 tensor ret(2, def_dim_2, 0.0); 00560 Cosseratstresstensor s( 0.0); 00561 Cosseratstresstensor t( 0.0); 00562 tensor I2("I", 2, def_dim_2); 00563 00564 // double EPS = pow(d_macheps(),(1./3.)); 00565 double J2D = this->Jinvariant2(); 00566 double q = this->q_deviatoric(); 00567 double theta = this->theta(); 00568 00569 //out while ( theta >= 2.0*PI ) 00570 //out theta = theta - 2.0*PI; // if bigger than full cycle 00571 //out while ( theta >= 4.0*PI/3.0 ) 00572 //out theta = theta - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00573 //out while ( theta >= 2.0*PI/3.0 ) 00574 //out theta = theta - 2.0*PI/3.0; // if bigger than two third of half cycle 00575 //out while ( theta >= PI/3.0 ) 00576 //out theta = 2.0*PI/3.0 - theta; // if bigger than one third of half cycle 00577 //out 00578 //out if ( theta < 0.0001 ) 00579 //out { 00580 //out ::printf("theta = %.10e in CAMPotentialSurface::dthetaoverds(stress)\n", 00581 //out theta); 00582 //out//>><<>><<>><< return ret; 00583 //out } 00584 00585 double c3t = cos(3.0*theta); 00586 double s3t = sin(3.0*theta); 00587 00588 double tempS = (3.0/2.0)*(c3t/(q*q*s3t)); 00589 double tempT = (9.0/2.0)*(1.0/(q*q*q*s3t)); 00590 00591 s = this->deviator(); 00592 t = s("qk")*s("kp") - I2*(J2D*(2.0/3.0)); 00593 00594 s.null_indices(); 00595 t.null_indices(); 00596 00597 ret = s*tempS - t*tempT; 00598 ret.null_indices(); 00599 return ret; 00600 } 00601 00602 //############################################################################# 00608 //############################################################################# 00609 tensor Cosseratstresstensor::d2poverds2( void ) const 00610 { 00611 tensor ret(4, def_dim_4, 0.0); 00612 return ret; 00614 } 00615 00616 00617 //############################################################################# 00618 tensor Cosseratstresstensor::d2qoverds2( void ) const 00619 { 00620 // Cosseratstresstensor stress = EPS->getStress(); 00621 00622 00623 //first the return tensor in order not to fragment the memory 00624 tensor ret( 4, def_dim_4, 0.0); // second derivative of q over 00625 // d sigma_pq d sigma_mn 00626 //setting up some constants 00627 00628 tensor I2("I", 2, def_dim_2); // To check out this three fourth-order 00629 tensor I_pqmn("I", 4, def_dim_4); // isotropic tensor please see 00630 I_pqmn = I2("pq")*I2("mn"); // W.Michael Lai, David Rubin, 00631 I_pqmn.null_indices(); 00632 tensor I_pmqn("I", 4, def_dim_4); // Erhard Krempl 00633 I_pmqn = I_pqmn.transpose0110(); // " Introduction to Continuum Mechanics" 00634 // QA808.2 ; ISBN 0-08-022699-X 00635 00636 double q_dev = this->q_deviatoric(); 00637 00638 Cosseratstresstensor s = this->deviator(); 00639 s.null_indices(); 00640 00641 tensor iso1( 4, def_dim_4, 0.0); //this will be d_pm*d_nq-d_pq*d_nm*(1/3) 00642 iso1 = I_pmqn - I_pqmn*(1.0/3.0); 00643 00644 double tempiso1 = (3.0/2.0)*(1/q_dev); 00645 double tempss = (9.0/4.0)*(1.0/( q_dev * q_dev * q_dev )); 00646 00647 ret = iso1*tempiso1 - (s("pq")*s("mn"))*tempss; 00648 ret.null_indices(); 00649 return ret; 00650 } 00651 00652 00653 //############################################################################# 00654 tensor Cosseratstresstensor::d2thetaoverds2( void ) const 00655 { 00656 tensor ret( 4, def_dim_4, 0.0); 00657 00658 tensor I2("I", 2, def_dim_2); 00659 tensor I_pqmn("I", 4, def_dim_4); // To check out this three fourth-order 00660 I_pqmn = I2("pq")*I2("mn"); // isotropic tensor please see 00661 I_pqmn.null_indices(); 00662 tensor I_pmqn("I", 4, def_dim_4); // W.Michael Lai, David Rubin, 00663 I_pmqn = I_pqmn.transpose0110(); // Erhard Krempl 00664 //no tensor I_pnqm("I", 4, def_dim_4); // " Introduction to Continuum Mechanics" 00665 //no I_pnqm = I_pqmn.transpose0111(); // QA808.2 ; ISBN 0-08-022699-X 00666 00667 double J2D = this->Jinvariant2(); 00668 00669 // double EPS = pow(d_macheps(),(1./3.)); 00670 00671 tensor s( 2, def_dim_2, 0.0); 00672 tensor t( 2, def_dim_2, 0.0); 00673 s = this->deviator(); 00674 t = s("qk")*s("kp") - I2*(J2D*(2.0/3.0)); 00675 //s.print("s"," \n\n tensor s \n\n"); 00676 //t.print("t"," \n\n tensor t \n\n"); 00677 s.null_indices(); 00678 t.null_indices(); 00679 00680 tensor p( 4, def_dim_4, 0.0); //this will be d_mp*d_nq - d_pq*d_nm*(1/3) 00681 tensor w( 4, def_dim_4, 0.0); 00682 00683 double theta = this->theta(); 00684 //out while ( theta >= 2.0*PI ) 00685 //out theta = theta - 2.0*PI; // if bigger than full cycle 00686 //out while ( theta >= 4.0*PI/3.0 ) 00687 //out theta = theta - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00688 //out while ( theta >= 2.0*PI/3.0 ) 00689 //out theta = theta - 2.0*PI/3.0; // if bigger than two third of half cycle 00690 //out while ( theta >= PI/3.0 ) 00691 //out theta = 2.0*PI/3.0 - theta; // if bigger than one third of half cycle 00692 //out 00693 //out if ( theta < 0.0001 ) 00694 //out { 00695 //out ::printf("\n\ntheta = %.10e in CAMPotentialSurface::d2thetaoverds2(stress)\n", 00696 //out theta); 00697 //out//>><<>><<>><< return ret; 00698 //out } 00699 //out 00700 00701 double q_dev = this->q_deviatoric(); 00702 00703 //setting up some constants 00704 //...... double sqrt3 = sqrt(3.0); 00705 double c3t = cos(3*theta); 00706 double s3t = sin(3*theta); 00707 double s3t3 = s3t*s3t*s3t; 00708 double q3 = q_dev * q_dev * q_dev; 00709 double q4 = q_dev * q_dev * q_dev * q_dev; 00710 double q5 = q_dev * q_dev * q_dev * q_dev * q_dev; 00711 double q6 = q_dev * q_dev * q_dev * q_dev * q_dev * q_dev; 00712 00713 00714 double tempss = - (9.0/2.0)*(c3t)/(q4*s3t) - (27.0/4.0)*(c3t/(s3t3*q4)); 00715 00716 double tempst = (81.0/4.0)*(1.0)/(s3t3*q5); 00717 00718 double tempts = (81.0/4.0)*(1.0/(s3t*q5)) + (81.0/4.0)*(c3t*c3t)/(s3t3*q5); 00719 00720 double temptt = - (243.0/4.0)*(c3t/(s3t3*q6)); 00721 00722 double tempp = +(3.0/2.0)*(c3t/(s3t*q_dev*q_dev)); 00723 00724 double tempw = -(9.0/2.0)*(1.0/(s3t*q3)); 00725 00726 00727 //)))))))) 00728 //)))))))) double tempss = 0.0; 00729 //)))))))) 00730 //)))))))) double tempst = 0.0; 00731 //)))))))) 00732 //)))))))) double tempts = 0.0; 00733 //)))))))) 00734 //)))))))) double temptt = 0.0; 00735 //)))))))) 00736 //)))))))) double tempp = 0.0; 00737 //)))))))) 00738 //)))))))) double tempw = 0.0; 00739 //)))))))) 00740 00741 00742 // fourth order tensors in the final equation for second derivative 00743 // of theta over ( d \sigma_{pq} d \sigma_{mn} ) 00744 // BE CAREFULL order is PQ MN 00745 00746 tensor s_pq_d_mn( 4, def_dim_4, 0.0); 00747 tensor s_pn_d_mq( 4, def_dim_4, 0.0); 00748 //... tensor s_pm_d_nq( 4, def_dim_4, 0.0); 00749 00750 tensor d_pq_s_mn( 4, def_dim_4, 0.0); 00751 tensor d_pn_s_mq( 4, def_dim_4, 0.0); 00752 //... tensor d_pm_s_nq( 4, def_dim_4, 0.0); 00753 00754 p = I_pmqn - I_pqmn*(1.0/3.0); 00755 00756 // privremena stampa 00757 //.......................................................................... 00758 //.......p.print("p"," tensor p = I_pmqn - I_pqmn*(1.0/3.0)"); 00759 00760 00761 s_pq_d_mn = s("pq")*I2("mn"); 00762 s_pq_d_mn.null_indices(); 00763 s_pn_d_mq = s_pq_d_mn.transpose0101(); 00764 //... s_pm_d_nq = s_pn_d_mq.transpose0110(); 00765 00766 d_pq_s_mn = I2("pq")*s("mn"); 00767 d_pq_s_mn.null_indices(); 00768 d_pn_s_mq = d_pq_s_mn.transpose0101(); 00769 //... d_pm_s_nq = d_pn_s_mq.transpose0110(); 00770 00773 //s_pq_d_mn.print("sd"," tensor s_pq_d_mn = s(\"pq\")*I2(\"mn\") "); 00774 //s_pn_d_mq.print("sd"," tensor s_pn_d_mq = s_pq_d_mn.transpose0101()"); 00776 //d_pq_s_mn.print("ds"," tensor d_pq_s_mn = I2(\"pq\")*s(\"mn\") "); 00777 //d_pn_s_mq.print("ds"," tensor d_pn_s_mq = d_pq_s_mn.transpose0101()"); 00780 00781 w = s_pn_d_mq + d_pn_s_mq - s_pq_d_mn*(2.0/3.0) - d_pq_s_mn*(2.0/3.0); 00782 00783 //.....// symmetric w 00784 //..... w = ( s_pn_d_mq + s_pm_d_nq ) * 0.5 + 00785 //..... ( d_pn_s_mq + d_pm_s_nq ) * 0.5 00786 //..... - s_pq_d_mn*(2.0/3.0) - d_pq_s_mn*(2.0/3.0); 00787 //..... 00788 00791 //w.print("w","w=s_pn_d_mq+d_pn_s_mq-s_pq_d_mn*(2.0/3.0)-d_pq_s_mn*(2.0/3.0)"); 00793 00796 //tensor ss = s("pq")*s("mn"); 00797 //tensor st = s("pq")*t("mn"); 00798 //tensor ts = t("pq")*s("mn"); 00799 //tensor tt = t("pq")*t("mn"); 00800 //ss.print("ss","\n tensor ss "); 00801 //st.print("st","\n tensor st "); 00802 //ts.print("ts","\n tensor ts "); 00803 //tt.print("tt","\n tensor tt "); 00805 00806 // finally 00807 ret = ( s("pq")*s("mn") * tempss 00808 + s("pq")*t("mn") * tempst 00809 + t("pq")*s("mn") * tempts 00810 + t("pq")*t("mn") * temptt 00811 + p * tempp 00812 + w * tempw ); 00813 00814 ret.null_indices(); 00815 return ret; 00816 } 00817 00818 00819 //--//############################################################################## 00820 //--// trying to find intersection point 00821 //--// according to M. Crisfield's book 00822 //--// " Non-linear Finite Element Analysis of Solids and Structures " 00823 //--// chapter 6.6.1 page 168. 00824 //--Cosseratstresstensor Cosseratstresstensor::yield_surface_cross(Cosseratstresstensor & end_stress, 00825 //-- Material_Model & YC) 00826 //--{ 00827 //--// bounds 00828 //-- double x1 = 0.0; 00829 //-- double x2 = 1.0; 00830 //--// accuracy 00831 //-- double const TOL = 1.0e-9; 00832 //-- double a = zbrentstress( *this, end_stress, YC, x1, x2, TOL ); 00833 //--// ::printf("\n****\n a = %lf \n****\n",a); 00834 //-- Cosseratstresstensor delta_stress = end_stress - *this; 00835 //-- Cosseratstresstensor intersection_stress = *this + delta_stress * a; 00836 //--//*** intersection_this->reportshort("FINAL Intersection stress\n"); 00837 //-- 00838 //-- return intersection_stress; 00839 //-- 00840 //--} 00841 00842 //############################################################################## 00843 Cosseratstresstensor Cosseratstresstensor::pqtheta2stress( double p, double q, double theta) 00844 { 00845 Cosseratstresstensor ret; 00846 // ret.report("ret"); 00847 // double sqrtJ2D = q/3.0; 00848 // double 2osqrt3 = 2.0/sqrt(3.0); 00849 double temp = q*TWOOVERTHREE; 00850 00851 // while ( theta >= 2.0*PI ) 00852 // theta = theta - 2.0*PI; // if bigger than full cycle 00853 // while ( theta >= 4.0*PI/3.0 ) 00854 // theta = theta - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00855 // while ( theta >= 2.0*PI/3.0 ) 00856 // theta = theta - 2.0*PI/3.0; // if bigger than two third of half cycle 00857 // while ( theta >= PI/3.0 ) 00858 // theta = 2.0*PI/3.0 - theta; // if bigger than one third of half cycle 00859 00860 double ct = cos( theta ); 00861 double ctm = cos( theta - TWOOVERTHREE*PI ); 00862 double ctp = cos( theta + TWOOVERTHREE*PI ); 00863 00864 ret.val(1,1) = temp*ct - p; 00865 ret.val(2,2) = temp*ctm - p; 00866 ret.val(3,3) = temp*ctp - p; 00867 00868 // ret.report("ret"); 00869 return ret; 00870 00871 } 00872 00873 //############################################################################## 00874 void Cosseratstresstensor::report(char * msg) const 00875 { 00876 ::printf("**************** stress tensor report ****************\n"); 00877 if ( msg ) ::printf("%s",msg); 00878 00879 this->print("st","Cosseratstresstensor st"); 00880 00881 ::fprintf(stdout,"I1 = %.8e ; I2 = %.8e ; I3 = %.8e ;\n", 00882 Iinvariant1(),Iinvariant2(),Iinvariant3()); 00883 00884 fprintf(stdout,"st_trace = %.8e, mean pressure p = %.8e\n ", 00885 trace(), trace()/3.0); 00886 00887 BJtensor I2("I", 2, def_dim_2); 00888 00889 Cosseratstresstensor st_vol = I2 * trace() * (1./3.); 00890 st_vol.print("st_v","tensor st_vol (volumetric part of the st tensor)"); 00891 00892 Cosseratstresstensor st_dev = this->deviator(); // - st_vol; 00893 st_dev.print("st_d","tensor st_dev (stress deviator)"); 00894 00895 ::fprintf(stdout,"J1 = %.8e ; J2 = %.8e ; J3 = %.8e ;\n", 00896 Jinvariant1(),Jinvariant2(),Jinvariant3()); 00897 //... Cosseratstresstensor Jinv2 = st_dev("ij")*st_dev("ji")*0.5; 00898 00899 Cosseratstresstensor st_principal = principal(); 00900 st_principal.print("st_p","principal stress tensor"); 00901 00902 00903 ::fprintf(stdout,"sig_oct = %.8e , tau_oct = %.8e\n", 00904 sigma_octahedral(), tau_octahedral()); 00905 00906 00907 ::printf("ksi = %.8e , ro = %.8e , theta = %.8e\n", 00908 ksi(), ro(), theta()); 00909 00910 ::printf("p=%.12e , q=%.12e , theta=%.12e*PI\n", 00911 p_hydrostatic(), q_deviatoric(), thetaPI()); 00912 00913 if ( msg ) ::printf("%s",msg); 00914 ::printf("############ end of stress tensor report ############\n"); 00915 } 00916 00917 00918 //############################################################################## 00919 void Cosseratstresstensor::reportshort(char * msg) const 00920 { 00921 //.. ::printf("\n** short stress tensor report ****************\n"); 00922 // if ( msg ) ::printf("** %s",msg); 00923 00924 this->print("st"," "); 00925 // ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00926 // ksi(), ro(), theta()); 00927 // ::printf(" p=%.8e q=%.8e theta=%.8e\n", 00928 // p_hydrostatic(), q_deviatoric(), theta() ); 00929 } 00930 00931 //############################################################################## 00932 void Cosseratstresstensor::reportshortpqtheta(char * msg) const 00933 { 00934 //.. ::printf("\n** short stress tensor report ****************\n"); 00935 if ( msg ) ::printf("%s",msg); 00936 00937 //.. this->print("st","Cosseratstresstensor st"); 00938 //.. ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00939 //.. ksi(), ro(), theta()); 00940 00941 double p = p_hydrostatic(); 00942 double q = q_deviatoric(); 00943 double t = theta(); 00944 ::printf(" p= %+.6e q= %+.6e theta= %+.6e \n", p, q, t); 00945 } 00946 //############################################################################## 00947 void Cosseratstresstensor::reportSHORTpqtheta(char * msg) const 00948 { 00949 //.. ::printf("\n** short stress tensor report ****************\n"); 00950 if ( msg ) ::printf("%s",msg); 00951 00952 //.. this->print("st","Cosseratstresstensor st"); 00953 //.. ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00954 //.. ksi(), ro(), theta()); 00955 00956 double p = p_hydrostatic(); 00957 double q = q_deviatoric(); 00958 double t = theta(); 00959 ::printf(" p= %+.6e q= %+.6e theta= %+.6e ", p, q, t); 00960 } 00961 //############################################################################## 00962 void Cosseratstresstensor::reportSHORTs1s2s3(char * msg) const 00963 { 00964 //.. ::printf("\n** short stress tensor report ****************\n"); 00965 if ( msg ) ::printf("%s",msg); 00966 00967 //.. this->print("st","Cosseratstresstensor st"); 00968 //.. ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00969 //.. ksi(), ro(), theta()); 00970 00971 this->print("st",""); 00972 } 00973 //############################################################################## 00974 void Cosseratstresstensor::reportKLOTpqtheta(char * msg) const 00975 { 00976 //.. ::printf("\n** short stress tensor report ****************\n"); 00977 if ( msg ) ::printf("%s",msg); 00978 00979 //.. this->print("st","Cosseratstresstensor st"); 00980 //.. ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00981 //.. ksi(), ro(), theta()); 00982 00983 double p = p_hydrostatic(); 00984 double q = q_deviatoric(); 00985 double t = theta(); 00986 ::printf(" %+.6e %+.6e %+.6e ", p, q, t); 00987 } 00988 //############################################################################## 00989 void Cosseratstresstensor::reportshortI1J2J3(char * msg) const 00990 { 00991 //.. ::printf("\n** short stress tensor report ****************\n"); 00992 if ( msg ) ::printf("%s",msg); 00993 00994 //.. this->print("st","Cosseratstresstensor st"); 00995 //.. ::printf("** ksi = %.8e , ro = %.8e , theta = %.8e\n", 00996 //.. ksi(), ro(), theta()); 00997 00998 double I1 = Iinvariant1(); 00999 double J2 = Jinvariant2(); 01000 double J3 = Jinvariant3(); 01001 ::printf(" I1= %.6e J2= %.6e J3= %.6e \n", I1, J2, J3); 01002 } 01003 //############################################################################## 01004 void Cosseratstresstensor::reportAnim(void) const 01005 { 01006 ::printf("Anim p= %f ; q= %f ; theta= %f \n", 01007 p_hydrostatic(), q_deviatoric(), theta()); 01008 } 01009 //############################################################################## 01010 void Cosseratstresstensor::reportTensor(char * msg) const 01011 { 01012 //.. ::printf("\n** short stress tensor report ****************\n"); 01013 if ( msg ) ::printf("%s",msg); 01014 01015 ::fprintf(stdout," %+.6e %+.6e %+.6e %+.6e %+.6e %+.6e \n", 01016 this->cval(1,1), 01017 this->cval(1,2), 01018 this->cval(1,3), 01019 this->cval(2,2), 01020 this->cval(2,3), 01021 this->cval(3,3)); 01022 } 01023 01024 01025 //############################################################################## 01026 OPS_Stream& operator<< (OPS_Stream& os, const Cosseratstresstensor & rhs) 01027 //ostream& operator<< (ostream& os, const Cosseratstresstensor & rhs) 01028 { 01029 //if ( msg ) ::printf("%s",msg); 01030 01031 // os.setf( ios::showpos | ios::scientific); 01032 01033 os.precision(4); 01034 os.width(10); 01035 01036 //os << endlnn; 01037 //os << rhs.cval(1,1) << " "; 01038 //os.width(10); 01039 //os << rhs.cval(1,2) << " "; 01040 //os.width(10); 01041 //os << rhs.cval(1,3) << endlnn; 01042 //os << rhs.cval(2,1) << " "; 01043 //os.width(10); 01044 //os << rhs.cval(2,2) << " "; 01045 //os.width(10); 01046 //os << rhs.cval(2,3) << endlnn; 01047 // 01048 //os << rhs.cval(3,1) << " "; 01049 //os.width(10); 01050 //os << rhs.cval(3,2) << " "; 01051 //os.width(10); 01052 //os << rhs.cval(3,3) << endlnn; 01053 01054 os.width(10); 01055 //os << "p = " << rhs.p_hydrostatic(); 01056 os << " "<< rhs.p_hydrostatic(); 01057 os.width(10); 01058 //os << " q = " << rhs.q_deviatoric(); 01059 os << " " << rhs.q_deviatoric(); 01060 os.width(10); 01061 //os << " theta = " << rhs.theta(); 01062 01063 return os; 01064 } 01065 01066 01067 01068 01069 #endif 01070 |