求教fish 编写遗传算法程序

ysj6174

如题 qq510384742

gao9941

这个问题很高深

yangxin546

我只会用FISH编写粒子群算法- -

zpp080425

fish编写遗传算法在FLAC当中不可能实现，FLAC不能自己调用自己，所以放手去找其他出路吧，这条路行不通。

chenpf09

好高深，，，，，，，，，               

jl66682008

const char *UserCamClayModel::Run(unsigned uDim,State *ps) {

if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserCamClayModel");

static double dPav,dQav,dEvav,dEvdPav;

static int iPind;

/* --- plasticity indicator:
---*/

/*
store 'now' info. as 'past' and turn 'now' info off ---*/

if (ps->mState & mShearNow) ps->mState |= mShearPast;

ps->mState &= ~mShearNow;

/* --- initialize stacks --- */

if (!ps->bySubZone) {

iPind = 0;

dPav
= 0.0;

dQav
= 0.0;

dEvav
= 0.0;

dEvdPav = 0.0;

}

/* --- trial elastic stresses试验弹性应力 --- */
double dA1 = dBulk + dC4D3 * dShear;
//体积模量K + 4.0 / 3.0*剪切模量G， dA1 =K+
double dA2 = dBulk - dC2D3 * dShear;
// dA2 =K-
double dG2 = 2.0 * dShear;
// dG2 =2

/*--由主应变分量增量得主应力增量 --*/


ps->stnS.d11 += ps->stnE.d11 * dA1 + (ps->stnE.d22 + ps->stnE.d33) * dA2

ps->stnS.d22 += ps->stnE.d22 * dA1 + (ps->stnE.d11 + ps->stnE.d33) * dA2;

ps->stnS.d33 += ps->stnE.d33 * dA1 + (ps->stnE.d11 + ps->stnE.d22) * dA2;
/* 因为在 中 */
ps->stnS.d12 += ps->stnE.d12 * dG2;

ps->stnS.d13 += ps->stnE.d13 * dG2;

ps->stnS.d23 += ps->stnE.d23 * dG2;

/* --- mean pressure平均压力dPVal
=
--- */

double dPVal = - (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) * dC1D3;

/* --- deviatoric stresses偏应力 --- */

double dDS11 = ps->stnS.d11 + dPVal;

double dDS22 = ps->stnS.d22 + dPVal;

double dDS33 = ps->stnS.d33 + dPVal;

double dDS12 = ps->stnS.d12;//

double dDS13 = ps->stnS.d13;

double dDS23 = ps->stnS.d23;
/*second deviatoric stress invariant第二偏应力不变量 - */
Double dJ2 =0.5*(dDS11*dDS11+dDS22*dDS22+dDS33*dDS33) + ps->stnS.d12* ps->stnS.d12 + ps->stnS.d13 * ps->stnS.d13 + ps->stnS.d23 * ps->stnS.d23;
/*----广义剪应力 ----*/

double dQVal = sqrt(3.*dJ2);

/* --- Cam-clay yield criterion剑桥屈服准则--- */

double dM2 = dMM*dMM;
//dMM是摩擦常数M
dM2 =M2

double dFS = dQVal*dQVal + dM2 * dPVal * (dPVal - dMPC); /* ----- 先期固结压力在这里是硬化参数Ha.即Ha=p c , f = ------*/

double dPNew = 0.0;

double dQNew = 0.0;

double dEVVal = 0.0;

double dEvdPVal = 0.0;

/* --- detect small pc --- */
if (dMPC < dMP1*1.e-3) {
/*先期固结压力<基准压力 *10-3 的情况下：*/

iPind = 1;

ps->mState |= mShearNow;

ps->stnS.d11 = 0.0;

// Sx=0

ps->stnS.d22 = 0.0;

ps->stnS.d33 = 0.0;

ps->stnS.d12 = 0.0;

ps->stnS.d13 = 0.0;

ps->stnS.d23 = 0.0;

dPNew
= dPVal;

dQNew
= dQVal;

dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=

dEvdPVal = dEVVal;

} else {

if (dFS > 0.0) {

//如果屈服函数大于零，表示剪切破坏

/* --- shear failure剪切破坏 --- */

iPind = 1;

ps->mState |= mShearNow;

double dSa
= 6.0 * dShear * dQVal;

//dSa=

double dSc
= dM2 * (2.0 * dPVal - dMPC);

//dSc=

double dSb
= dBulk * dSc;

//dSb =弹性体积模量
K*ca

double dBa
= dSa*dSa + dM2 * dSb*dSb;
//dBa=
double dBb = dSa * dQVal + 0.5 * dSb * dSc; //dBb =6qG*q+0.5* K* M2*(2p- pc) *M2*(2p- pc)=

double dBc
= dFS;
// dBc
=

double dAlam
= 0.0;

double dAlam1 = 0.0;

if (dBa != 0.0) {

double dBoa
= dBb / dBa;

//dBoa=
double dVal =dBoa*dBoa - dBc / dBa;
//dVal= 因为a ,b ,c 满足式 ，
if (dVal < 0.0) return("Cam-clay:
Yield envelope cannot be reached");//未达到屈服轨迹线

double dVal1
= sqrt(dVal);
//dVal1=


dAlam
= dBoa + dVal1;
// dAlam= +

dAlam1 = dBoa - dVal1;

//dAlam1= -

} else {

if (dBb != 0.0) {
dAlam
= 0.5 * dBc / dBb;
// dAlam=0.5* /
dAlam1 = dAlam;

} else return("Cam-clay:
Yield envelope cannot be reached");

}

if (fabs(dAlam1) < fabs(dAlam)) dAlam = dAlam1;

if (dAlam < 0.0) dAlam = 0.0;

dQNew
= dQVal - dSa * dAlam;
//

dPNew
= dPVal - dSb * dAlam;
//

/* --- update stresses 更新应力--- */

if (dQVal == 0.0) {


//如果 为0

double dAdd = dPNew - dPVal;
//

ps->stnS.d11 -= dAdd;

ps->stnS.d22 -= dAdd;

ps->stnS.d33 -= dAdd;

} else { //如果 不为0,新的应力分量： ,其中，偏应力分量：
ps->stnS.d11 = (dDS11 / dQVal) * dQNew - dPNew;
//




ps->stnS.d22 = (dDS22 / dQVal) * dQNew - dPNew;

ps->stnS.d33 = (dDS33 / dQVal) * dQNew - dPNew;

ps->stnS.d12 = (dDS12 / dQVal) * dQNew;

//

ps->stnS.d23 = (dDS23 / dQVal) * dQNew;

ps->stnS.d13 = (dDS13 / dQVal) * dQNew;

}

dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=

dEvdPVal = dAlam * dSc;

} else {

/* --- no failure --- */

dPNew
= dPVal;

dQNew
= dQVal;

dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=

dEvdPVal = 0.0;

}

}

/* update zone stresses and strains更新单元的应力应变--- */

double dVol = ps->dSubZoneVolume;

dPav
+= dPNew * dVol;
//在前面有/dPav,dQav,dEvav,dEvdPav;且设置初始值为0

dQav
+= dQNew * dVol;

dEvav
+= dEVVal * dVol;

dEvdPav += dEvdPVal * dVol;

/* --- the last zone has been processed, update parameters: 更新变量--- */

if (ps->bySubZone==ps->byTotSubZones-1) {

dVol = 1.0 / ps->dZoneVolume;

if (ps->byOverlay==2) dVol *= 0.5;

dPav
= dPav * dVol;

dQav
= dQav * dVol;

dEvav
= dEvav * dVol;

dEvdPav = dEvdPav * dVol;
dMP
= dPav;

//当前的平均有效应力

dMQ
= dQav;

//当前的平均偏应力

dEV
+= dEvav;
//累计总体积应变+=

dEV_P += dEvdPav;
//累计塑性体积应变

/* --- specific dVolume --- */

double dVold
= dMV;
//初始比容v0

dMV
= dVold * (1.0 + dEvav);

if (dMV < 0.0) return("Cam-clay:
Specific volume is negative");
//初始比容v0

/* ---

bulk modulus --- */

dBulk
= dMV * dPav / dKappa; //弹性体积模量K =初始比容v0*弹性膨胀线的斜率
if (dBulk > dBulkB) return("Cam-clay: Bulk modulus bound must be increased for stability");//弹性体积模量K
/* --- cap pressure硬化—软化准则：屈服曲线的形状取决于固结压力pc的值，此压力是塑性体积变化的函数，且随比容而变化。对于新的 , , ,可以通过 平面内固结线和膨胀线的交点得出，其表达式如下：--- */

if (iPind == 1) {

if (dPav > dMPC * 0.4) {
//先期固结压力

double dVal = dMV + dKappa * log(dPav / dMP1); /*初始比容v0+
弹性膨胀线的斜率
/基准压力 , 式中   


dMPC = dMP1 * exp((dMV_L - dVal)/(dLambda-dKappa)); /*---先期固结压力=基准压力 ╳e （正常固结线的比容 - 弹性膨胀线的斜率 - ---*/

} else {
dMPC=dMPC*(1.0+dEvdPav*dVold/(dLambda-dKappa));//先期固结压力=弹性膨胀线斜率

}

if (dMPC <= 0.0) { //先期固结压力小于等于

dMPC
= 0.0;
//那设置先期固结压力0，因为土不能抗拉

dBulk = 0.0;
//弹性体积模量K

}

}

/* --- shear modulus --- */

if (dPoisson != 0.0) {
//泊松比不等于0
dShear=1.5*dBulk*(1.0-2.0*dPoisson)/(1.0+dPoisson); //剪切模量, 弹性体积模量,泊松比
/* ,弹性体积模量 随比容和平均应力的函数变化的*/

} else {

double dMaxg = 1.5 * dBulk;
//弹性体积模量K

double dMing = 0.0;

dShear = dShear < dMaxg ? dShear : dMaxg; //剪切模量

dShear = dShear > dMing ? dShear : dMing;

}

}

return(0);
}

jl66682008

不好意思，粘错了