copy_con.cpp

本文介绍了一个使用C++实现的动态数组类,包括简单构造函数、析构函数及拷贝构造函数的实现。该类通过new操作符分配内存,并在拷贝构造函数中实现了深拷贝,确保了数组数据的正确复制。

  name="google_ads_frame" marginwidth="0" marginheight="0" src="http://pagead2.googlesyndication.com/pagead/ads?client=ca-pub-5572165936844014&dt=1194442938015&lmt=1194190197&format=336x280_as&output=html&correlator=1194442937843&url=file%3A%2F%2F%2FC%3A%2FDocuments%2520and%2520Settings%2Flhh1%2F%E6%A1%8C%E9%9D%A2%2FCLanguage.htm&color_bg=FFFFFF&color_text=000000&color_link=000000&color_url=FFFFFF&color_border=FFFFFF&ad_type=text&ga_vid=583001034.1194442938&ga_sid=1194442938&ga_hid=1942779085&flash=9&u_h=768&u_w=1024&u_ah=740&u_aw=1024&u_cd=32&u_tz=480&u_java=true" frameborder="0" width="336" scrolling="no" height="280" allowtransparency="allowtransparency"> #include <iostream.h>
#include <stdlib.h>

class array
{
   int *p;
   int size;
 public:
   array(int sz) {           // simple constructor
      p = new int[sz];
      if(!p) exit(1);
      size = sz;
    }
   ~array() {delete [] p;}    // destructor
   array(const array &object);// copy constructor
   void put(int i, int j){
      if(i>=0 && i<size)
        p[i] = j;
    }
   int get(int i) {return p[i];}
 };

array::array(const array &object)
 {
   int lcl_i;

   p = new int[object.size];
   if (!p)
      exit(1);
   for(lcl_i=0; lcl_i < object.size; lcl_i++)
      p[lcl_i] = object.p[lcl_i];
 }

void main(void)
 {
   array num(10);
   int lcl_i;

   for (lcl_i=0; lcl_i<10; lcl_i++)
      num.put(lcl_i, lcl_i);
   for (lcl_i=9; lcl_i>=0; lcl_i--)
      cout << num.get(lcl_i);
   cout << endl;

   //  Create another array using the copy constructor
   array x=num;
   for (lcl_i=0; lcl_i<10; lcl_i++)
      cout << x.get(lcl_i);
 }

 

 

#include <opencv2/opencv.hpp> #include <cmath> // 假设 solvePoisson 函数已存在 cv::Mat solvePoisson(const cv::Mat &rhs); cv::Mat unwrap_TIE(const cv::Mat &phase_wrap) { // 确保输入是单通道浮点矩阵 CV_Assert(phase_wrap.type() == CV_32FC1); int rows = phase_wrap.rows; int cols = phase_wrap.cols; // 步骤1: 计算 psi = exp(1i * phase_wrap) cv::Mat cos_phase, sin_phase; cv::cos(phase_wrap, cos_phase); // 实部: cos(phase) cv::sin(phase_wrap, sin_phase); // 虚部: sin(phase) // 步骤2: 计算x方向差分并包裹到[-π, π] cv::Mat diff_x = cv::Mat::zeros(rows, cols - 1, CV_32F); for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols - 1; ++j) { // 计算相邻像素的相位差: angle(conj(psi(i,j)) * psi(i,j+1)) float a_real = cos_phase.at<float>(i, j); float a_imag = sin_phase.at<float>(i, j); float b_real = cos_phase.at<float>(i, j + 1); float b_imag = sin_phase.at<float>(i, j + 1); // 计算 conj(a) * b float real_part = a_real * b_real + a_imag * b_imag; float imag_part = a_real * b_imag - a_imag * b_real; // 计算相位差并包裹到[-π, π] diff_x.at<float>(i, j) = std::atan2(imag_part, real_part); } } // 构建edx矩阵 (左右各填充一列0) cv::Mat edx = cv::Mat::zeros(rows, cols + 1, CV_32F); diff_x.copyTo(edx(cv::Rect(1, 0, cols - 1, rows))); // 中间部分填充相位差 // 步骤3: 计算y方向差分并包裹到[-π, π] cv::Mat diff_y = cv::Mat::zeros(rows - 1, cols, CV_32F); for (int i = 0; i < rows - 1; ++i) { for (int j = 0; j < cols; ++j) { // 计算相邻像素的相位差: angle(conj(psi(i,j)) * psi(i+1,j)) float a_real = cos_phase.at<float>(i, j); float a_imag = sin_phase.at<float>(i, j); float b_real = cos_phase.at<float>(i + 1, j); float b_imag = sin_phase.at<float>(i + 1, j); // 计算 conj(a) * b float real_part = a_real * b_real + a_imag * b_imag; float imag_part = a_real * b_imag - a_imag * b_real; // 计算相位差并包裹到[-π, π] diff_y.at<float>(i, j) = std::atan2(imag_part, real_part); } } // 构建edy矩阵 (上下各填充一行0) cv::Mat edy = cv::Mat::zeros(rows + 1, cols, CV_32F); diff_y.copyTo(edy(cv::Rect(0, 1, cols, rows - 1))); // 中间部分填充相位差 // 步骤4: 计算拉普拉斯算子 (二阶差分) cv::Mat dxx = edx(cv::Rect(1, 0, cols, rows)) - edx(cv::Rect(0, 0, cols, rows)); cv::Mat dyy = edy(cv::Rect(0, 1, cols, rows)) - edy(cv::Rect(0, 0, cols, rows)); cv::Mat lap = dxx + dyy; // 拉普拉斯算子 // 步骤5: 计算 rho = imag(conj(psi) .* lap) cv::Mat rho = -sin_phase.mul(lap); // 等价于 imag(conj(psi)*lap // 步骤6: 求解泊松方程 cv::Mat phase_unwrap = solvePoisson(rho); return phase_unwrap; } cv没有成员con,sin,编译con,sin。
07-18
#include "modelDeemo.h" #include "version.txt" // excerpt-export-start int __stdcall DllMain(void *,unsigned, void *) { return 1; } extern "C" EXPORT_TAG const char *getName() { #ifdef MODELDEBUG return "cmodelDeemod"; #else return "cmodelDeemo"; #endif } extern "C" EXPORT_TAG unsigned getMajorVersion() { return MAJOR_VERSION; } extern "C" EXPORT_TAG unsigned getMinorVersion() { return UPDATE_VERSION; } extern "C" EXPORT_TAG void *createInstance() { models::ModelDeemo *m = new models::ModelDeemo(); return (void *)m; } // excerpt-export-end namespace models { // excerpt-con-start ModelDeemo::ModelDeemo() : bulk_(0.0), Mshear_(0.0), Mviscosity1_(0.0), Kshear1_(0.0), Kviscosity1_(0.0), Kshear2_(0.0), Kviscosity2_(0.0), Mviscosity2_(0.0), A_(0.0) { //kelvinStrain1 Mekd_[0] = 0.0; Mekd_[1] = 0.0; Mekd_[2] = 0.0; Mekd_[3] = 0.0; Mekd_[4] = 0.0; Mekd_[5] = 0.0; //kelvinStrain2 Mekd_[6] = 0.0; Mekd_[7] = 0.0; Mekd_[8] = 0.0; Mekd_[9] = 0.0; Mekd_[10] = 0.0; Mekd_[11] = 0.0; } // excerpt-con-end String ModelDeemo::getName() const { #ifdef MODELDEBUG return L"Deemo-debug"; #else return L"Deemo"; #endif } String ModelDeemo::getFullName() const { #ifdef MODELDEBUG return L"Deemo Debug"; #else return L"Deemo"; #endif } UInt ModelDeemo::getMinorVersion() const { return UPDATE_VERSION; } String ModelDeemo::getProperties() const { //获取参数 return L"bulk_, Kshear1_, Mshear_, Kviscosity1_, Mviscosity1_,"//burgers的参数 L"Kshear2_, Kviscosity2_, Mviscosity2_, A_," //塑性部分的参数 L"strain-kelvin-xx,strain-kelvin-yy,strain-kelvin-zz,strain-kelvin-xy,strain-kelvin-xz,strain-kelvin-yz,"//黏弹性体应变 L"cohesion,friction,dilation,tension,strain-shear-plastic,strain-tensile-plastic -strain-tension-plastic";//MC准则 } String ModelDeemo::getStates() const { return L"shear-n,tension-n,shear-p,tension-p"; } Variant ModelDeemo::getProperty(UInt index) const { switch (index) { case 1: return bulk_; case 2: return Kshear1_; case 3: return Mshear_; case 4: return Kviscosity1_; case 5: return Mviscosity1_; case 6: return Kshear2_; case 7: return Kviscosity2_; case 8: return Mviscosity2_; case 9: return A_;//模型参数 case 10: return Mekd_[0]; case 11: return Mekd_[1]; case 12: return Mekd_[2]; case 13: return Mekd_[3]; case 14: return Mekd_[4]; case 15: return Mekd_[5]; case 16: return Mekd_[6]; case 17: return Mekd_[7]; case 18: return Mekd_[8]; case 19: return Mekd_[9]; case 20: return Mekd_[10]; case 21: return Mekd_[11]; case 22: return cohesion_; case 23: return friction_; case 24: return dilation_; case 25: return tension_; case 26: return sHP_;//剪切应变 case 27: return tHP_;//张拉应变 } return(0.0); } void ModelDeemo::setProperty(UInt index,const Variant &p,UInt restoreVersion) { ConstitutiveModel::setProperty(index, p, restoreVersion); switch (index) { case 1: bulk_ = p.toDouble(); break; case 2: Kshear1_ = p.toDouble(); break; case 3: Mshear_ = p.toDouble(); break; case 4: Kviscosity1_ = p.toDouble(); break; case 5: Mviscosity1_ = p.toDouble(); break; case 6: Kshear2_ = p.toDouble(); break; case 7: Kviscosity2_ = p.toDouble(); break; case 8: Mviscosity2_ = p.toDouble(); break; case 9: A_ = p.toDouble(); break; case 10: Mekd_[0] = p.toDouble(); break; case 11: Mekd_[1] = p.toDouble(); break; case 12: Mekd_[2] = p.toDouble(); break; case 13: Mekd_[3] = p.toDouble(); break; case 14: Mekd_[4] = p.toDouble(); break; case 15: Mekd_[5] = p.toDouble(); break; case 16: Mekd_[6] = p.toDouble(); break; case 17: Mekd_[7] = p.toDouble(); break; case 18: Mekd_[8] = p.toDouble(); break; case 19: Mekd_[9] = p.toDouble(); break; case 20: Mekd_[10] = p.toDouble(); break; case 21: Mekd_[11] = p.toDouble(); break; case 22: cohesion_ = p.toDouble(); break; case 23: friction_ = p.toDouble(); break; case 24: dilation_ = p.toDouble(); break; case 25: tension_ = p.toDouble(); break; case 26: sHP_ = p.toDouble(); break; case 27: tHP_ = p.toDouble(); break; } } /*bool ModelDeemo::isPropertyAdvanced(UInt i) const { if (i <= 4) return ModelDeemo::isPropertyAdvanced(i); else if (i==9) return true; return false; }*/ //下面这段cvisc中是vm,这里是mm void ModelDeemo::copy(const ConstitutiveModel *m) { const ModelDeemo *mm = dynamic_cast<const ModelDeemo *>(m); if (!mm) throw std::runtime_error("Internal error: constitutive model dynamic cast failed."); // ConstitutiveModel::copy(m); // bulk_ = mm->bulk_; Kshear1_ = mm->Kshear1_; Mshear_ = mm->Mshear_; Kviscosity1_ = mm->Kviscosity1_; Mviscosity1_ = mm->Mviscosity1_; Kshear2_ = mm->Kshear2_; Kviscosity2_ = mm->Kviscosity2_; Mviscosity2_ = mm->Mviscosity2_; A_ = mm->A_; Mekd_[0] = mm->Mekd_[0]; Mekd_[1] = mm->Mekd_[1]; Mekd_[2] = mm->Mekd_[2]; Mekd_[3] = mm->Mekd_[3]; Mekd_[4] = mm->Mekd_[4]; Mekd_[5] = mm->Mekd_[5]; Mekd_[6] = mm->Mekd_[6]; Mekd_[7] = mm->Mekd_[7]; Mekd_[8] = mm->Mekd_[8]; Mekd_[9] = mm->Mekd_[9]; Mekd_[10] = mm->Mekd_[10]; Mekd_[11] = mm->Mekd_[11]; cohesion_ = mm->cohesion_; friction_ = mm->friction_; dilation_ = mm->dilation_; tension_ = mm->tension_; sHP_ = mm->sHP_; tHP_ = mm->tHP_; } // excerpt-run-start参数初始化,计算模型中间参数值; void ModelDeemo::initialize(UByte dim, State *s) { ConstitutiveModel::initialize(dim, s); //参数非零化 if (Mshear_ <= 0.0) Mshear_ = 1e-20; if (Kshear1_ <= 0.0) Kshear1_ = 0.0; if (Kshear2_ <= 0.0) Kshear2_ = 0.0; if (Kviscosity1_ <= 0.0) { Kviscosity1_ = 0.0; Kshear1_ = 0.0; } if (Kviscosity2_ <= 0.0) { Kviscosity2_ = 0.0; Kshear2_ = 0.0; } if (Mviscosity1_ <= 0.0) Mviscosity1_ = 0.0; if (Mviscosity2_ <= 0.0) Mviscosity2_ = 0.0; if (A_ <= 0.0) A_ = 0.0; // 3. 初始化Mohr-Coulomb准则相关参数(参考CVisc模型) Double rsin = std::sin(friction_ * degrad); // degrad是角度转弧度的系数 nph_ = (1.0 + rsin) / (1.0 - rsin); // 塑性硬化参数 csn_ = 2.0 * cohesion_ * sqrt(nph_); // 修正后的内聚力 // 限制抗拉强度不超过内聚力和摩擦角决定的顶点值 if (friction_ > 0.0) { Double apex = cohesion_ / std::tan(friction_ * degrad); tension_ = std::min(tension_, apex); } // 剪胀角相关参数 rsin = std::sin(dilation_ * degrad); nps_ = (1.0 + rsin) / (1.0 - rsin); // 剪胀参数 // 计算应力比参数 rc_ = std::sqrt(1.0 + nph_ * nph_); } static const UInt Dqs = 12; static const UInt Dqt = 13; void ModelDeemo::run(UByte dim, State *s) { ConstitutiveModel::run(dim, s); // excerpt-state-start if (s->state_ & shear_now) s->state_ |= shear_past; s->state_ &= ~shear_now; if (s->state_ & tension_now) s->state_ |= tension_past; s->state_ &= ~tension_now; // excerpt-state-end UInt iPlas = 0;//未屈服(黏弹性) // dEkd values now stored in s->working_[] array (necessary for thread safety) Double tempk1 = 0, tempk2 = 0, tempm1 = 0, tempm2 = 0.0; Double dCrtdel = (s->isCreep() ? s->getTimeStep() : 0.0); //sum += dCrtdel / 10.0;///加的东西 Double dSubZoneVolume = s->getSubZoneVolume(); if (!s->sub_zone_) { s->working_[0] = 0.0; s->working_[1] = 0.0; s->working_[2] = 0.0; s->working_[3] = 0.0; s->working_[4] = 0.0; s->working_[5] = 0.0; s->working_[6] = 0.0; s->working_[7] = 0.0; s->working_[8] = 0.0; s->working_[9] = 0.0; s->working_[10] = 0.0; s->working_[11] = 0.0; s->working_[Dqs] = 0.0; s->working_[Dqt] = 0.0; } if (Kviscosity1_ <= 0.0) tempk1 = 0.0; else tempk1 = 1.0 / Kviscosity1_; // if (Mviscosity1_ <= 0.0) tempm1 = 0.0; else tempm1 = 1.0 / Mviscosity1_; if (Kviscosity2_ <= 0.0) tempk2 = 0.0; else tempk2 = 1.0 / Kviscosity2_; if (Mviscosity2_ <= 0.0) tempm2 = 0.0; else tempm2 = 1.0 / Mviscosity2_; Double temp1 = 0.5 * Kshear1_ * dCrtdel * tempk1; Double a1_con = 1.0 + temp1;//A1 Double b1_con = 1.0 - temp1;//B1 Double ba1 = b1_con / a1_con; Double bac1 = ba1 - 1.0;//(B1/A1-1) Double temp2 = 0.5 * Kshear2_ * dCrtdel * tempk2; Double a2_con = 1.0 + temp2;//A2 Double b2_con = 1.0 - temp2;//B2 Double ba2 = b2_con / a2_con; Double bac2 = ba2 - 1.0;//(B2/A2-1) Double tempb1 = (tempm1 + tempk1 / a1_con ) * dCrtdel * 0.25; Double tempa = 0.5 / Mshear_; Double x_con1 = tempa + tempb1;//a1 Double y_con1 = tempa - tempb1;//b1 Double c1dxc1 = 1.0 / x_con1;//1/a1 Double tempb2 = (tempm1 + tempk1 / a1_con+ tempk2 / a2_con) * dCrtdel * 0.25; Double x_con2 = tempa + tempb2;//a2 Double y_con2 = tempa - tempb2;//b2 Double c1dxc2 = 1.0 / x_con2;//1/a2 Double z1_con = dCrtdel * tempk1 / (4.0 * a1_con); Double z2_con = dCrtdel * tempk2 / (4.0 * a2_con); //Double z_con = z1_con + z2_con;// //;--- partition strains --- Double dev = s->stnE_.s11() + s->stnE_.s22() + s->stnE_.s33(); Double dev3 = d1d3 * dev; Double de11d = s->stnE_.s11() - dev3; Double de22d = s->stnE_.s22() - dev3; Double de33d = s->stnE_.s33() - dev3;//应变计算 //;--- partition stresses--- Double s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33()); Double s11d = s->stnS_.s11() - s0; Double s22d = s->stnS_.s22() - s0; Double s33d = s->stnS_.s33() - s0;//应力计算 //;--- remember old stresses --- Double s11old = s11d; Double s22old = s22d; Double s33old = s33d; Double s12old = s->stnS_.s12(); Double s13old = s->stnS_.s13(); Double s23old = s->stnS_.s23(); //;--- new trial deviator stresses assuming viscoelastic increments --- s11d = (de11d + s11d * y_con1 - Mekd_[0] * bac1) * c1dxc1; s22d = (de22d + s22d * y_con1 - Mekd_[1] * bac1) * c1dxc1; s33d = (de33d + s33d * y_con1 - Mekd_[2] * bac1) * c1dxc1; Double s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con1 - Mekd_[3] * bac1) * c1dxc1; Double s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con1 - Mekd_[4] * bac1) * c1dxc1; Double s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con1 - Mekd_[5] * bac1) * c1dxc1; //;--- new trial isotropic stress assuming elastic increment --- s0 += bulk_ * dev; Double s11i = s11d + s0; Double s22i = s22d + s0; Double s33i = s33d + s0; //; --- trial stresses --- s->stnS_.rs11() = s11i; s->stnS_.rs22() = s22i; s->stnS_.rs33() = s33i; s->stnS_.rs12() = s12i; s->stnS_.rs13() = s13i; s->stnS_.rs23() = s23i; //先标定了没有进入塑性的部分 //下面是在做塑性状态的判断 if (canFail()) { SymTensorInfo info; DVect3 prin = s->stnS_.getEigenInfo(&info); Double e1_ = bulk_ + d2d3 * c1dxc2;//α1 Double e2_ = bulk_ - d1d3 * c1dxc2;//α2 Double ra = e1_ - nps_ * e2_; //Double e21_ = e2_ / e1_; Double rb = e2_ * (1.0 - nps_); Double rd = e2_ - nps_ * e1_;//应力更新公式的λ的系数 //定义塑性部分相关参数 Double G3 = Kshear2_; // 开尔文体弹性模量 (对应Kshear2_) Double A = A_; // 非线性黏壶参数 Double fs = -prin.x() + nph_ * prin.z() - csn_; Double ftz = prin.z() - tension_; Double fty = prin.y() - tension_; Double ftx = prin.x() - tension_;//MC破坏准则 Double fsd = fs / rc_;//(归一化但是还没搞清楚是啥) Double lambda_s_1 = fs / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel)); //λs1 Double lambda_s_2 = fs * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); //λs2 Double lambda_s = lambda_s_1 + lambda_s_2;//λs Double beta_s_xt = ftx / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βx Double gamma_s_xt = ftx * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γx Double lambda_s_xt = beta_s_xt + gamma_s_xt;//λsxt Double beta_s_yt = fty / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βy Double gamma_s_yt = fty * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γy Double lambda_s_yt = beta_s_yt + gamma_s_yt;//λsyt Double beta_s_zt = ftz / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βz Double gamma_s_zt = ftz * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γz Double lambda_s_zt = beta_s_zt + gamma_s_zt;//λszt /*if (iPlas == 0) { // 塑性未激活:仅用第一个 Kelvin 元件(原逻辑) s11d = (de11d + s11d * y_con1 - Mekd_[0] * bac1) * c1dxc1; s22d = (de22d + s22d * y_con1 - Mekd_[1] * bac1) * c1dxc1; s33d = (de33d + s33d * y_con1 - Mekd_[2] * bac1) * c1dxc1; s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con1 - Mekd_[3] * bac1) * c1dxc1; s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con1 - Mekd_[4] * bac1) * c1dxc1; s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con1 - Mekd_[5] * bac1) * c1dxc1; } else { //Double s11_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s11d)) / (2 * m1_); //Double s22_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s22d)) / (2 * m1_); //Double s33_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s33d)) / (2 * m1_); //Double s12_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s12())) / (2 *m1_); //Double s13_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s13())) / (2 *m1_); //Double s23_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s23())) / (2 *m1_); //非线性牛顿体 Double s11_2 = s11d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); Double s22_2 = s22d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); Double s33_2 = s33d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); Double s12_2 = s->stnS_.s12() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); Double s13_2 = s->stnS_.s13() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); Double s23_2 = s->stnS_.s23() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); // 塑性已激活:使用两个 Kelvin 元件和一个非线性牛顿体 s11d = (de11d + s11d * y_con2 - s11_2 - Mekd_[0] * bac1 - Mekd_[6] * bac2) * c1dxc2; s22d = (de22d + s22d * y_con2 - s22_2 - Mekd_[1] * bac1 - Mekd_[7] * bac2) * c1dxc2; s33d = (de33d + s33d * y_con2 - s33_2 - Mekd_[2] * bac1 - Mekd_[8] * bac2) * c1dxc2; s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con2 - s12_2 - Mekd_[3] * bac1 - Mekd_[9] * bac2) * c1dxc2; s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con2 - s13_2 - Mekd_[4] * bac1 - Mekd_[10] * bac2) * c1dxc2; s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con2 - s23_2 - Mekd_[5] * bac1 - Mekd_[11] * bac2) * c1dxc2; }//应力更新 */ if (fsd > 0.0 && fsd >= ftz) {//剪切塑性修正 iPlas = 1; s->state_ |= shear_now; prin.rx() += lambda_s * ra; prin.ry() += lambda_s * rb; prin.rz() += lambda_s * rd;//已检查 } else if (ftz > 0.0 && ftz >= fsd) { s->state_ |= tension_now;//拉伸塑性修正 if (ftx > 0.0) { iPlas = 4; prin.rx() -= lambda_s_xt * e1_ + (lambda_s_yt + lambda_s_zt) * e2_; prin.ry() -= lambda_s_yt * e1_ + (lambda_s_xt + lambda_s_zt) * e2_; prin.rz() -= lambda_s_zt * e1_ + (lambda_s_xt + lambda_s_yt) * e2_;//三个主应力失效 } else if (fty > 0.0) { iPlas = 3; prin.rx() -= (lambda_s_yt + lambda_s_zt) * e2_; prin.ry() -= (lambda_s_zt * e2_ + lambda_s_yt * e1_); prin.rz() -= (lambda_s_zt * e1_ + lambda_s_yt * e2_);//两个主应力失效 } else { iPlas = 2; //Double tco = ftz * e21_;//应该是用不上了 prin.rx() -= lambda_s_zt * e2_ ; prin.ry() -= lambda_s_zt * e2_ ; prin.rz() -= lambda_s_zt * e1_ ;//一个主应力失效 } } if (friction_ > 0.0) { Double apex = cohesion_ / tan(friction_*degrad); if (prin.x() >= apex || prin.y() >= apex || prin.z() >= apex) {//如果任意方向上的应力超过这个阈值,进入拉伸状态 iPlas = 4; s->state_ |= tension_now; prin.rx() -= lambda_s_xt * e1_ + (lambda_s_yt + lambda_s_zt) * e2_; prin.ry() -= lambda_s_yt * e1_ + (lambda_s_xt + lambda_s_zt) * e2_; prin.rz() -= lambda_s_zt * e1_ + (lambda_s_xt + lambda_s_yt) * e2_; } } if (iPlas) { s->stnS_ = info.resolve(prin);//用当前的prin应力来解析并更新应力状态 if (iPlas == 1) { Double dDe1p = -lambda_s;//1方向的塑性应变增量 Double dDe3p = -1 * lambda_s * nps_;//3方向的塑性应变增量 Double dDepa = d1d3 * (dDe1p + dDe3p); dDe1p -= dDepa; dDe3p -= dDepa; s->working_[Dqs] += sqrt(0.5 * (dDe1p*dDe1p + dDepa * dDepa + dDe3p * dDe3p)) * dSubZoneVolume; } if (iPlas == 2) { Double dAux = lambda_s_zt; s->working_[Dqt] += dAux * dSubZoneVolume; } if (iPlas == 3) { Double dAux = lambda_s_zt + lambda_s_yt; s->working_[Dqt] += dAux * dSubZoneVolume; } if (iPlas == 4) { Double dAux = lambda_s_zt + lambda_s_yt + lambda_s_xt; s->working_[Dqt] += dAux * dSubZoneVolume; } } } //--- sub-zone contribution to Kelvin1-strains --- s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33()); s->working_[0] += (Mekd_[0] * ba1 + (s->stnS_.s11() - s0 + s11old) * z1_con) * s->getSubZoneVolume(); s->working_[1] += (Mekd_[1] * ba1 + (s->stnS_.s22() - s0 + s22old) * z1_con) * s->getSubZoneVolume(); s->working_[2] += (Mekd_[2] * ba1 + (s->stnS_.s33() - s0 + s33old) * z1_con) * s->getSubZoneVolume(); s->working_[3] += (Mekd_[3] * ba1 + (s->stnS_.s12() + s12old) * z1_con) * s->getSubZoneVolume(); s->working_[4] += (Mekd_[4] * ba1 + (s->stnS_.s13() + s13old) * z1_con) * s->getSubZoneVolume(); s->working_[5] += (Mekd_[5] * ba1 + (s->stnS_.s23() + s23old) * z1_con) * s->getSubZoneVolume(); //--- sub-zone contribution to Kelvin2-strains --- s->working_[6] += (Mekd_[6] * ba2 + (s->stnS_.s11() - s0 + s11old) * z2_con) * s->getSubZoneVolume(); s->working_[7] += (Mekd_[7] * ba2 + (s->stnS_.s22() - s0 + s22old) * z2_con) * s->getSubZoneVolume(); s->working_[8] += (Mekd_[8] * ba2 + (s->stnS_.s33() - s0 + s33old) * z2_con) * s->getSubZoneVolume(); s->working_[9] += (Mekd_[9] * ba2 + (s->stnS_.s12() + s12old) * z2_con) * s->getSubZoneVolume(); s->working_[10] += (Mekd_[10] * ba2 + (s->stnS_.s13() + s13old) * z2_con) * s->getSubZoneVolume(); s->working_[11] += (Mekd_[11] * ba2 + (s->stnS_.s23() + s23old) * z2_con) * s->getSubZoneVolume(); //--- update stored Kelvin-strains and plastic strain --- if (s->sub_zone_ == s->total_sub_zones_ - 1) {//这行代码检查是否在最后一个子区域内,如果是,则执行后续的更新操作 Double Aux = 1. / s->getZoneVolume();//Aux变量是区域体积的倒数 if (s->overlay_ == 2) Aux *= 0.5;//如果overlay_等于2,则Aux乘以0.5. Mekd_[0] = s->working_[0] * Aux; Mekd_[1] = s->working_[1] * Aux; Mekd_[2] = s->working_[2] * Aux; Mekd_[3] = s->working_[3] * Aux; Mekd_[4] = s->working_[4] * Aux; Mekd_[5] = s->working_[5] * Aux; Mekd_[6] = s->working_[6] * Aux; Mekd_[7] = s->working_[7] * Aux; Mekd_[8] = s->working_[8] * Aux; Mekd_[9] = s->working_[9] * Aux; Mekd_[10] = s->working_[10] * Aux; Mekd_[11] = s->working_[11] * Aux; if (canFail()) { sHP_ += s->working_[Dqs] * Aux; tHP_ += s->working_[Dqt] * Aux; } } // excerpt-run-end } } // namespace models // EOF具体到我这个本构你看一下应该怎么改呢,是不是前面的mekd6-11也应该删掉,void ModelCVisc::run(UByte d,State *s) { ConstitutiveModel::run(d,s); if (s->state_ & shear_now) s->state_ |= shear_past; s->state_ &= ~shear_now; if (s->state_ & tension_now) s->state_ |= tension_past; s->state_ &= ~tension_now; UInt iPlas = 0; Double tempk=0, tempm=0; Double dCrtdel = (s->isCreep() ? s->getTimeStep() : 0.0); Double dSubZoneVolume = s->getSubZoneVolume(); if (!s->sub_zone_) { s->working_[0] = 0.0; s->working_[1] = 0.0; s->working_[2] = 0.0; s->working_[3] = 0.0; s->working_[4] = 0.0; s->working_[5] = 0.0; s->working_[Dqs] =0.0; s->working_[Dqt] =0.0; } if(Kviscosity_ <= 0.0) tempk = 0.0; else tempk = 1.0 / Kviscosity_ ; // if (Mviscosity_ <= 0.0) tempm = 0.0; else tempm = 1.0 / Mviscosity_ ; Double temp = 0.5 * Kshear_ * dCrtdel * tempk; Double a_con = 1.0 + temp; Double b_con = 1.0 - temp; Double ba = b_con / a_con; Double bal = ba - 1.0; temp = (tempm + tempk / a_con) * dCrtdel * 0.25 ; Double temp1 = 0.5 / Mshear_; Double x_con = temp1 + temp; Double y_con = temp1 - temp; Double z_con = dCrtdel * tempk / (4.0 * a_con); Double c1dxc = 1.0 / x_con; //;--- partition strains --- Double dev = s->stnE_.s11() + s->stnE_.s22() + s->stnE_.s33() ; Double dev3 = d1d3 * dev; Double de11d = s->stnE_.s11() - dev3; Double de22d = s->stnE_.s22() - dev3; Double de33d = s->stnE_.s33() - dev3; //;--- partition stresses--- Double s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33()); Double s11d = s->stnS_.s11() - s0; Double s22d = s->stnS_.s22() - s0; Double s33d = s->stnS_.s33() - s0; //;--- remember old stresses --- Double s11old = s11d; Double s22old = s22d; Double s33old = s33d; Double s12old = s->stnS_.s12(); Double s13old = s->stnS_.s13(); Double s23old = s->stnS_.s23(); //;--- new trial deviator stresses assuming viscoelastic increments --- s11d = (de11d + s11d * y_con - Mekd_[0] * bal) * c1dxc ; s22d = (de22d + s22d * y_con - Mekd_[1]* bal) * c1dxc ; s33d = (de33d + s33d * y_con - Mekd_[2]* bal) * c1dxc ; Double s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - Mekd_[3] * bal) * c1dxc ; Double s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - Mekd_[4] * bal) * c1dxc ; Double s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - Mekd_[5] * bal) * c1dxc ; //;--- new trial isotropic stress assuming elastic increment --- s0 += bulk_ * dev; Double s11i = s11d + s0; Double s22i = s22d + s0; Double s33i = s33d + s0; //; --- trial stresses --- s->stnS_.rs11() = s11i; s->stnS_.rs22() = s22i; s->stnS_.rs33() = s33i; s->stnS_.rs12() = s12i; s->stnS_.rs13() = s13i; s->stnS_.rs23() = s23i; 前面的部分是不是要改成这样,这个是前面的Burgers模型的,我本构模型的前面部分是和他一样的,后面在进行我的代码里面的,应力修正这样是对的吗
10-30
#include "modelcvisc.h" #include "base/src/version.h" namespace models { ModelCVisc::ModelCVisc() : ModelBurger(), cohesion_(0.0),friction_(0.0),dilation_(0.0),tension_(0.0), sHP_(0.0),tHP_(0.0),nph_(0.0),csn_(0.0),rc_(0.0),nps_(0.0) { } UInt ModelCVisc::getMinorVersion() const{ return version::update; } String ModelCVisc::getProperties(void) const { return ModelBurger::getProperties() + L",cohesion,friction,dilation,tension,strain-shear-plastic,strain-tensile-plastic -strain-tension-plastic"; } String ModelCVisc::getStates(void) const { return L"shear-n,tension-n,shear-p,tension-p"; } Variant ModelCVisc::getProperty(UInt index) const { if (index <= 11) return ModelBurger::getProperty(index); else { switch (index) { case 12: return cohesion_; case 13: return friction_; case 14: return dilation_; case 15: return tension_; case 16: return sHP_; case 17: return tHP_; } } return(0.0); } void ModelCVisc::setProperty(UInt index, const Variant &p,UInt restoreVersion) { if (index <= 11) ModelBurger::setProperty(index,p,restoreVersion); else { switch (index) { case 12: cohesion_ = p.toDouble(); break; case 13: friction_ = p.toDouble(); break; case 14: dilation_ = p.toDouble(); break; case 15: tension_ = p.toDouble(); break; case 16: sHP_ = p.toDouble(); break; case 17: tHP_ = p.toDouble(); break; } } } void ModelCVisc::copy(const ConstitutiveModel *m) { const ModelCVisc *vm = dynamic_cast<const ModelCVisc *>(m); if (!vm) throw std::runtime_error("Internal error: constitutive model dynamic cast failed."); ModelBurger::copy(m); cohesion_ = vm->cohesion_; friction_ = vm->friction_; dilation_ = vm->dilation_; tension_ = vm->tension_; sHP_ = vm->sHP_; tHP_ = vm->tHP_; } void ModelCVisc::initialize(UByte d,State *s) { ModelBurger::initialize(d, s); Double rsin = std::sin(friction_ * degrad); nph_ = (1.0 + rsin) / (1.0 - rsin); csn_ = 2.0 * cohesion_ * sqrt(nph_); if (friction_) { Double apex = cohesion_ / std::tan(friction_ * degrad); tension_ = std::min(tension_,apex); } rsin = std::sin(dilation_ * degrad); nps_ = (1.0 + rsin) / (1.0 - rsin); rc_ = std::sqrt(1.0 + nph_*nph_); } static const UInt Dqs = 6; static const UInt Dqt = 7; void ModelCVisc::run(UByte d,State *s) { ConstitutiveModel::run(d,s); if (s->state_ & shear_now) s->state_ |= shear_past; s->state_ &= ~shear_now; if (s->state_ & tension_now) s->state_ |= tension_past; s->state_ &= ~tension_now; UInt iPlas = 0; Double tempk=0, tempm=0; Double dCrtdel = (s->isCreep() ? s->getTimeStep() : 0.0); Double dSubZoneVolume = s->getSubZoneVolume(); if (!s->sub_zone_) { s->working_[0] = 0.0; s->working_[1] = 0.0; s->working_[2] = 0.0; s->working_[3] = 0.0; s->working_[4] = 0.0; s->working_[5] = 0.0; s->working_[Dqs] =0.0; s->working_[Dqt] =0.0; } if(Kviscosity_ <= 0.0) tempk = 0.0; else tempk = 1.0 / Kviscosity_ ; // if (Mviscosity_ <= 0.0) tempm = 0.0; else tempm = 1.0 / Mviscosity_ ; Double temp = 0.5 * Kshear_ * dCrtdel * tempk; Double a_con = 1.0 + temp; Double b_con = 1.0 - temp; Double ba = b_con / a_con; Double bal = ba - 1.0; temp = (tempm + tempk / a_con) * dCrtdel * 0.25 ; Double temp1 = 0.5 / Mshear_; Double x_con = temp1 + temp; Double y_con = temp1 - temp; Double z_con = dCrtdel * tempk / (4.0 * a_con); Double c1dxc = 1.0 / x_con; //;--- partition strains --- Double dev = s->stnE_.s11() + s->stnE_.s22() + s->stnE_.s33() ; Double dev3 = d1d3 * dev; Double de11d = s->stnE_.s11() - dev3; Double de22d = s->stnE_.s22() - dev3; Double de33d = s->stnE_.s33() - dev3; //;--- partition stresses--- Double s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33()); Double s11d = s->stnS_.s11() - s0; Double s22d = s->stnS_.s22() - s0; Double s33d = s->stnS_.s33() - s0; //;--- remember old stresses --- Double s11old = s11d; Double s22old = s22d; Double s33old = s33d; Double s12old = s->stnS_.s12(); Double s13old = s->stnS_.s13(); Double s23old = s->stnS_.s23(); //;--- new trial deviator stresses assuming viscoelastic increments --- s11d = (de11d + s11d * y_con - Mekd_[0] * bal) * c1dxc ; s22d = (de22d + s22d * y_con - Mekd_[1]* bal) * c1dxc ; s33d = (de33d + s33d * y_con - Mekd_[2]* bal) * c1dxc ; Double s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - Mekd_[3] * bal) * c1dxc ; Double s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - Mekd_[4] * bal) * c1dxc ; Double s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - Mekd_[5] * bal) * c1dxc ; //;--- new trial isotropic stress assuming elastic increment --- s0 += bulk_ * dev; Double s11i = s11d + s0; Double s22i = s22d + s0; Double s33i = s33d + s0; //; --- trial stresses --- s->stnS_.rs11() = s11i; s->stnS_.rs22() = s22i; s->stnS_.rs33() = s33i; s->stnS_.rs12() = s12i; s->stnS_.rs13() = s13i; s->stnS_.rs23() = s23i; if (canFail()) { SymTensorInfo info; DVect3 prin = s->stnS_.getEigenInfo(&info); Double e1_ = bulk_ + d2d3*c1dxc; Double e2_ = bulk_ - d1d3*c1dxc; Double ra = e1_ - nps_*e2_; Double rb = e2_ - nps_*e1_; Double rd = ra - rb*nph_; Double sc1_ = ra / rd; Double sc3_ = rb / rd; Double sc2_ = e2_*(1.0 - nps_)/rd; Double e21_ = e2_/e1_; Double sf1_ = 1.0/rd; Double sf3_ = - nps_/rd; Double fs = - prin.x() + nph_ * prin.z() - csn_; Double fsd = fs / rc_; Double ftz = prin.z() - tension_; Double fty = prin.y() - tension_; Double ftx = prin.x() - tension_; if (fsd > 0.0 && fsd >= ftz) { iPlas = 1; s->state_ |= shear_now; prin.rx() += fs * sc1_; prin.ry() += fs * sc2_; prin.rz() += fs * sc3_; } else if (ftz > 0.0 && ftz >= fsd) { s->state_ |= tension_now; if (ftx > 0.0) { iPlas = 4; prin.rx() = tension_; prin.ry() = tension_; prin.rz() = tension_; } else if (fty > 0.0) { iPlas = 3; prin.rx() -= (fty + ftz) * e2_/(e1_+e2_); prin.ry() = tension_; prin.rz() = tension_; } else { iPlas = 2; Double tco = ftz * e21_; prin.rx() -= tco; prin.ry() -= tco; prin.rz() = tension_; } } if (friction_ > 0.0) { Double apex = cohesion_ / tan(friction_*degrad); if (prin.x()>=apex || prin.y()>=apex || prin.z()>=apex) { iPlas = 4; s->state_ |= tension_now; prin.rx() = apex; prin.ry() = apex; prin.rz() = apex; } } if (iPlas) { s->stnS_ = info.resolve(prin); if (iPlas == 1) { Double dDe1p = fs * sf1_; Double dDe3p = fs * sf3_; Double dDepa = d1d3 * (dDe1p + dDe3p); dDe1p -= dDepa; dDe3p -= dDepa; s->working_[Dqs] += sqrt(0.5 * (dDe1p*dDe1p+dDepa*dDepa+dDe3p*dDe3p)) * dSubZoneVolume; } if (iPlas == 2) { Double dAux = ftz / e1_; s->working_[Dqt] += dAux * dSubZoneVolume; } if (iPlas == 3) { Double dAux = (ftz+fty) / (e1_+e2_); s->working_[Dqt] += dAux * dSubZoneVolume; } if (iPlas == 4) { Double dAux = (ftz+fty+ftx) / bulk_; s->working_[Dqt] += dAux * dSubZoneVolume; } } } //--- sub-zone contribution to Kelvin-strains --- s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33()); s->working_[0] += (Mekd_[0] * ba + (s->stnS_.s11() - s0 + s11old) * z_con) * s->getSubZoneVolume(); s->working_[1] += (Mekd_[1] * ba + (s->stnS_.s22() - s0 + s22old) * z_con) * s->getSubZoneVolume(); s->working_[2] += (Mekd_[2] * ba + (s->stnS_.s33() - s0 + s33old) * z_con) * s->getSubZoneVolume(); s->working_[3] += (Mekd_[3] * ba + (s->stnS_.s12() + s12old) * z_con) * s->getSubZoneVolume(); s->working_[4] += (Mekd_[4] * ba + (s->stnS_.s13() + s13old) * z_con) * s->getSubZoneVolume(); s->working_[5] += (Mekd_[5] * ba + (s->stnS_.s23() + s23old) * z_con) * s->getSubZoneVolume(); //--- update stored Kelvin-strains and plastic strain --- if (s->sub_zone_ == s->total_sub_zones_-1) { Double Aux = 1./s->getZoneVolume(); if (s->overlay_==2) Aux *= 0.5; Mekd_[0]= s->working_[0] * Aux; Mekd_[1]= s->working_[1] * Aux; Mekd_[2]= s->working_[2] * Aux; Mekd_[3]= s->working_[3] * Aux; Mekd_[4]= s->working_[4] * Aux; Mekd_[5]= s->working_[5] * Aux; if (canFail()) { sHP_ += s->working_[Dqs]*Aux; tHP_ += s->working_[Dqt]*Aux; } } } } // EOF这个是burgers-mohr的,看看他是这样实现的嘛
10-23
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值