OpenSees » SRC » material » nD » soil

MultiYieldSurface.cpp

Go to the documentation of this file.
00001 // $Revision: 1.8 $
00002 // $Date: 2003/02/14 23:01:30 $
00003 // $Source: /usr/local/cvs/OpenSees/SRC/material/nD/soil/MultiYieldSurface.cpp,v $
00004                                                                         
00005 // Written: ZHY
00006 // Created: August 2000
00007 //
00008 // MultiYieldSurface.cpp
00009 // ---------------------
00010 //
00011 
00012 #include <math.h>
00013 #include <stdlib.h>
00014 
00015 #include <MultiYieldSurface.h>
00016 
00017 
00018 // YieldSurface class methods
00019 MultiYieldSurface::MultiYieldSurface():
00020 theSize(0.0), theCenter(6), plastShearModulus(0.0)
00021 {
00022 
00023 }
00024 
00025 MultiYieldSurface::MultiYieldSurface(const Vector & theCenter_init, 
00026                                      double theSize_init, double plas_modul):
00027 theSize(theSize_init), theCenter(theCenter_init), plastShearModulus(plas_modul)
00028 {
00029 
00030 }
00031 
00032 MultiYieldSurface::~MultiYieldSurface()
00033 {
00034 
00035 }
00036 
00037 void MultiYieldSurface::setData(const Vector & theCenter_init, 
00038                                 double theSize_init, double plas_modul)
00039 {
00040   theSize = theSize_init;
00041   theCenter = theCenter_init;
00042   plastShearModulus = plas_modul;
00043 }
00044 
00045 void MultiYieldSurface::setCenter(const Vector & newCenter)
00046 {
00047   if (newCenter.Size() != 6) {
00048     opserr << "FATAL:MultiYieldSurface::setCenter(Vector &): vector size not equal 6" << endln;
00049     exit(-1);
00050   }
00051 
00052   theCenter = newCenter;
00053 }
00054 
00055 
00056 /**********************************************
00057 ostream & operator<< (ostream & os, const MultiYieldSurface & a)
00058 {
00059   os << "  theSize = " << a.theSize << endln 
00060      << "  theCenter = " << a.theCenter << endln
00061      << "  plastShearModulus = " << a.plastShearModulus << endln;
00062   
00063   return os;
00064 }
00065 
00066 
00067 istream & operator>> (istream & is, MultiYieldSurface & a)
00068 {
00069   is >> a.theSize >> a.theCenter >> a.plastShearModulus;
00070 
00071   return is;
00072 }
00073 *********************************************/
00074 
00075 double secondOrderEqn(double A, double B, double C, int i)
00076 {
00077   if(A == 0){
00078     opserr << "FATAL:second_order_eqn: A=0." << endln;
00079     if(i==0) opserr << " when finding reference point on outer surface." <<endln;
00080     else opserr << " when moving active surface." <<endln;
00081     exit(-1);   
00082   }
00083   if(C == 0) return 0;
00084   if(B == 0){
00085     if(C/A > 0){
00086       opserr << "FATAL:second_order_eqn: Complex roots.\n";
00087       exit(-1);
00088     } 
00089     return sqrt(-C/A);
00090   }
00091 
00092         double determ, val1, val2, val;
00093   determ = B*B - 4.*A*C; 
00094   if(determ < 0){
00095     opserr << "FATAL:second_order_eqn: Complex roots.\n";
00096     if(i==0) opserr << " when finding reference point on outer surface." <<endln;
00097     else opserr << " when moving active surface." <<endln;
00098     opserr << "B2=" << B*B << " 4AC=" << 4.*A*C <<endln; 
00099     exit(-1);
00100   }
00101   
00102   if (B > 0) val1 = (-B - sqrt(determ)) / (2.*A);
00103   else val1 = (-B + sqrt(determ)) / (2.*A);
00104   val2 = C / (A * val1);
00105 
00106   if (val1 < 0 && val2 < 0){
00107                 if (fabs(val1) < LOW_LIMIT) val1 = 0.;
00108                 else if (fabs(val2) < LOW_LIMIT) val2 = 0.;
00109         }
00110 
00111   if (val1 < 0 && val2 < 0){
00112     opserr << "FATAL:second_order_eqn: Negative roots.\n";
00113     if(i==0) opserr << " when finding reference point on outer surface." <<endln;
00114     else opserr << " when moving active surface." <<endln;
00115                 opserr << "A=" << A << " B=" << B << " C=" << C << " det=" << determ << 
00116                         " x1=" << val1 << " x2=" << val2 << endln;  
00117     exit(-1);   
00118   }
00119   
00120   if (val1 < 0) return  val2;
00121   else if (val2 < 0) return  val1;
00122   else{
00123     val = val1;
00124     if(val > val2)  val = val2;
00125                 return val;
00126   }
00127 }

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