C++ Matrix3::Exponential方法代码示例

// Increment small strain material deformation gradient using F(n) = (I + grad u)F(n-1)
void Elastic::HypoIncrementDeformation(MPMBase *mptr,Matrix3 du) const
{
	// get incremental deformation gradient
	const Matrix3 dF = du.Exponential(1);
	
	// current deformation gradient
	Matrix3 pF = mptr->GetDeformationGradientMatrix();
	
	// new deformation matrix
	const Matrix3 F = dF*pF;
    mptr->SetDeformationGradientMatrix(F);
}

// Entry point for large rotation
void IsoPlasticity::LRConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np,void *properties,ResidualStrains *res) const
{
	// current previous deformation gradient and stretch
	Matrix3 pFnm1 = mptr->GetDeformationGradientMatrix();
	
    // get incremental deformation gradient and decompose it
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
    Matrix3 dR;
    Matrix3 dV = dF.LeftDecompose(&dR,NULL);
	
	// decompose to get previous stretch
	Matrix3 Vnm1 = pFnm1.LeftDecompose(NULL,NULL);
	
	// get strain increments de = (dV-I) dR Vnm1 dRT
	dV(0,0) -= 1.;
	dV(1,1) -= 1.;
	dV(2,2) -= 1.;
	Matrix3 de = dV*Vnm1.RMRT(dR);
	
	// Update total deformation gradient
	Matrix3 pF = dF*pFnm1;
	mptr->SetDeformationGradientMatrix(pF);
	
    // Effective strain by deducting thermal strain (no shear thermal strain because isotropic)
	//  (note: using unreduced terms in CTE3 and CME3)
	double eres=CTE3*res->dT;
	if(DiffusionTask::active)
		eres+=CME3*res->dC;
	
	// Trial update assuming elastic response
	double delV;
    
    // 3D or 2D
	PlasticProperties *p = (PlasticProperties *)properties;
    if(np==THREED_MPM)
    {   delV = de.trace() - 3.*eres;
        LRPlasticityConstLaw(mptr,de(0,0),de(1,1),de(2,2),2*de(0,1),2.*de(0,2),2.*de(1,2),
                           delTime,np,delV,eres,p,res,&dR);
        return;
    }
	else if(np==PLANE_STRESS_MPM)
		delV = p->psRed*(de(0,0)+de(1,1)-2.*eres);
	else
		delV = de.trace() - 3.*eres;
    LRPlasticityConstLaw(mptr,de(0,0),de(1,1),2.*de(0,1),de(2,2),delTime,np,delV,eres,p,res,&dR);
}

/*  Given matrix of incremental deformation dF = exp(dt*grad v), increment particle strain,
        rotation, and LeftCauchy Green strain (latter is assumed to be stored in the particle's
        plastic strain tensor (which is accessed also with GetAltStrainTensor().
    New new F is dF.F, which is used to find new strain
    New B = dF.(Old B).dF^T
    Returns |dF|
*/
double HyperElastic::IncrementDeformation(MPMBase *mptr,Matrix3 du,Tensor *Btrial,int np) const
{
    // get incremental deformation gradient
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
	
	// current deformation gradient
	Matrix3 pF = mptr->GetDeformationGradientMatrix();
	
	// new deformation matrix
	const Matrix3 F = dF*pF;
    mptr->SetDeformationGradientMatrix(F);
	
    // increment Left Cauchy tensor B = F.F^T = dF.old B.dF^T
    // plain stress will need to update B.zz when known
    Matrix3 pBold = mptr->GetElasticLeftCauchyMatrix();
    
    // elements of dF.B
    Matrix3 dFoldB = dF*pBold;
    
    // return trial B (if provided) or store new B on the particle
    Tensor *pB = Btrial!=NULL ? Btrial : mptr->GetAltStrainTensor() ;
    pB->xx = dFoldB(0,0)*dF(0,0) + dFoldB(0,1)*dF(0,1) + dFoldB(0,2)*dF(0,2);
    pB->xy = dFoldB(0,0)*dF(1,0) + dFoldB(0,1)*dF(1,1) + dFoldB(0,2)*dF(1,2);
    
    pB->yy = dFoldB(1,0)*dF(1,0) + dFoldB(1,1)*dF(1,1) + dFoldB(1,2)*dF(1,2);
    
    pB->zz = dFoldB(2,0)*dF(2,0) + dFoldB(2,1)*dF(2,1) + dFoldB(2,2)*dF(2,2);
    
    if(np == THREED_MPM)
    {   pB->xz = dFoldB(0,0)*dF(2,0) + dFoldB(0,1)*dF(2,1) + dFoldB(0,2)*dF(2,2);
        pB->yz = dFoldB(1,0)*dF(2,0) + dFoldB(1,1)*dF(2,1) + dFoldB(1,2)*dF(2,2);
    }
    
    // return |dF|
    return dF.determinant();
}

// Apply Constitutive law, check np to know what type
void TaitLiquid::MPMConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np,void *properties,
									ResidualStrains *res,int historyOffset) const
{
#ifdef NO_SHEAR_MODEL
    // get incremental deformation gradient
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
    
    // decompose dF into dR and dU
    Matrix3 dR;
    Matrix3 dU = dF.RightDecompose(&dR,NULL);
	
	// current deformation gradient
    double detdF = dF.determinant();
	Matrix3 pF = mptr->GetDeformationGradientMatrix();
    Matrix3 F = dR*pF;
    if(np==THREED_MPM)
        F.Scale(pow(detdF,1./3.));
    else
        F.Scale2D(sqrt(detdF));
	
	// new deformation matrix with volume change onle
    mptr->SetDeformationGradientMatrix(F);

#else
#ifdef ELASTIC_B_MODEL
    // get incremental deformation gradient
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
	double detdF = dF.determinant();
	
	// current deformation gradient
	Matrix3 pF = mptr->GetDeformationGradientMatrix();
	
	// new deformation matrix
	const Matrix3 F = dF*pF;
    mptr->SetDeformationGradientMatrix(F);
#else

	// Update total deformation gradient, and calculate trial B
	double detdF = IncrementDeformation(mptr,du,NULL,np);
#endif
#endif
    
    // Get new J and save result on the particle
	double J = detdF * mptr->GetHistoryDble(J_History,historyOffset);
    mptr->SetHistoryDble(J_History,J,historyOffset);

#ifdef ELASTIC_B_MODEL
	// store pressure strain as elastic B
    Tensor *pB = mptr->GetAltStrainTensor() ;
    if(np==THREED_MPM || np==AXISYMMETRIC_MPM)
	{	double J23 = pow(J,2./3.);
		pB->xx = J23;
		pB->yy = J23;
		pB->zz = J23;
	}
	else
	{	pB->xx = J;
		pB->yy = J;
	}
#endif
    
    // account for residual stresses
	double dJres = GetIncrementalResJ(mptr,res);
	double Jres = dJres * mptr->GetHistoryDble(J_History+1,historyOffset);
	mptr->SetHistoryDble(J_History+1,Jres,historyOffset);
    double Jeff = J/Jres;

    // new Kirchhoff pressure (over rho0) from Tait equation
	double p0=mptr->GetPressure();
    double pressure = J*TAIT_C*Ksp*(exp((1.-Jeff)/TAIT_C)-1.);
    mptr->SetPressure(pressure);
    
	// incremental energy per unit mass - dilational part
    double avgP = 0.5*(p0+pressure);
    double delV = 1. - 1./detdF;
    double workEnergy = -avgP*delV;
    
	// incremental residual energy per unit mass
	double delVres = 1. - 1./dJres;
	double resEnergy = -avgP*delVres;
	
    // viscosity term = 2 eta (0.5(grad v) + 0.5*(grad V)^T - (1/3) tr(grad v) I) = 2 eta dev(grad v)
    // (i.e., deviatoric part of the symmetric strain tensor, 2 is for conversion to engineering shear strain)
	// simple shear rate = |2 dev(grad v)|
    Matrix3 shear;
    double c[3][3];
	double shearRate;
    c[0][0] = (2.*du(0,0)-du(1,1)-du(2,2))/3.;
    c[1][1] = (2.*du(1,1)-du(0,0)-du(2,2))/3.;
    c[2][2] = (2.*du(2,2)-du(0,0)-du(1,1))/3.;
    c[0][1] = 0.5*(du(0,1)+du(1,0));
    c[1][0] = c[0][1];
	shearRate = c[0][0]*c[0][0] + c[1][1]*c[1][1] + c[2][2]*c[2][2]
	+ 2.*c[0][1]*c[0][1];
    if(np==THREED_MPM)
    {   c[0][2] = 0.5*(du(0,2)+du(2,0));
        c[2][0] = c[0][2];
        c[1][2] = 0.5*(du(1,2)+du(2,1));
        c[2][1] = c[1][2];
//.........这里部分代码省略.........

/* Take increments in strain and calculate new Particle: strains, rotation strain,
	stresses, strain energy,
	du are (gradient rates X time increment) to give deformation gradient change
	For Axisymmetry: x->R, y->Z, z->theta, np==AXISYMMETRIC_MPM, otherwise dvzz=0
	This material tracks pressure and stores deviatoric stress only in particle stress tensor
 */
void Viscoelastic::MPMConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np,void *properties,ResidualStrains *res,int historyOffset) const
{
	// current previous deformation gradient and stretch
	Matrix3 pFnm1 = mptr->GetDeformationGradientMatrix();
	
    // get incremental deformation gradient and decompose it
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
    Matrix3 dR;
    Matrix3 dVstretch = dF.LeftDecompose(&dR,NULL);
	
	// decompose to get previous stretch
	Matrix3 Vnm1 = pFnm1.LeftDecompose(NULL,NULL);
	
	// get strain increments de = (dV-I) dR Vnma dRT
	dVstretch(0,0) -= 1.;
	dVstretch(1,1) -= 1.;
	dVstretch(2,2) -= 1.;
	Matrix3 Vrot = Vnm1.RMRT(dR);
	Matrix3 detot = dVstretch*Vrot;
	
	// Update total deformation gradient
	Matrix3 pF = dF*pFnm1;
	mptr->SetDeformationGradientMatrix(pF);
	
    // Effective strain by deducting thermal strain (no shear thermal strain because isotropic)
	double eres=CTE*res->dT;
	if(DiffusionTask::active)
		eres+=CME*res->dC;
	
	// update pressure
	double dTq0 = 0.,dispEnergy = 0.;
	double traceDe = detot.trace();
	double delV = traceDe - 3.*eres;
#ifdef USE_KIRCHOFF_STRESS
	// tracking J
	double detdF = dF.determinant();
	double **h =(double **)(mptr->GetHistoryPtr(0));
	double J = detdF*h[mptrHistory][MGJ_HISTORY];
	h[mptrHistory][MGJ_HISTORY] = J;
	UpdatePressure(mptr,delV,res,eres,detdF,J,delTime,dTq0,dispEnergy);
#else
	UpdatePressure(mptr,delV,res,eres,&dF,delTime,dTq0,dispEnergy);
#endif
	
	// deviatoric strains increment in de
	// Actually de is finding 2*(dev e) to avoid many multiplies by two
	Tensor de;
	double dV = traceDe/3.;
	de.xx = 2.*(detot(0,0) - dV);
	de.yy = 2.*(detot(1,1) - dV);
	de.zz = 2.*(detot(2,2) - dV);
	de.xy = 2.*detot(0,1);
	if(np==THREED_MPM)
	{	de.xz = 2.*detot(0,2);
		de.yz = 2.*detot(1,2);
	}
	
	// Find initial 2*e(t) (deviatoric strain) in ed
	Tensor ed;
	double thirdV = Vrot.trace()/3.;
	ed.xx = 2.*(Vrot(0,0)-thirdV);
	ed.yy = 2.*(Vrot(1,1)-thirdV);
	ed.zz = 2.*(Vrot(2,2)-thirdV);
	ed.xy = 2.*Vrot(0,1);
	if(np==THREED_MPM)
	{	ed.xz = 2.*Vrot(0,2);
		ed.yz = 2.*Vrot(1,2);
	}
	
	// increment particle deviatoric stresses - elastic part
	double dsig[6];
	dsig[XX] = Gered*de.xx;
	dsig[YY] = Gered*de.yy;
	dsig[ZZ] = Gered*de.zz;
	dsig[XY] = Gered*de.xy;
	if(np==THREED_MPM)
	{	dsig[XZ] = Gered*de.xz;
		dsig[YZ] = Gered*de.yz;
	}
	
	// get internal variable increments, update them, add to incremental stress, and get dissipated energy6
	Tensor dak;
    double **ak =(double **)(mptr->GetHistoryPtr(0));
	int k;
    for(k=0;k<ntaus;k++)
    {   double tmp = exp(-delTime/tauk[k]);
		double tmpm1 = tmp-1.;
		double tmpp1 = tmp+1.;
		double arg = 0.25*delTime/tauk[k];					// 0.25 because e's have factor of 2
		dak.xx = tmpm1*ak[XX_HISTORY][k] + arg*(tmpp1*ed.xx + de.xx);
		dak.yy = tmpm1*ak[YY_HISTORY][k] + arg*(tmpp1*ed.yy + de.yy);
		dak.xy = tmpm1*ak[XY_HISTORY][k] + arg*(tmpp1*ed.xy + de.xy);
		dak.zz = tmpm1*ak[ZZ_HISTORY][k] + arg*(tmpp1*ed.zz + de.zz);
		ak[XX_HISTORY][k] += dak.xx;
//.........这里部分代码省略.........

// Apply Constitutive law, check np to know what type
void TaitLiquid::MPMConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np,void *properties,ResidualStrains *res) const
{
#ifdef NO_SHEAR_MODEL
    // get incremental deformation gradient
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
    
    // decompose dF into dR and dU
    Matrix3 dR;
    Matrix3 dU = dF.RightDecompose(&dR,NULL);
	
	// current deformation gradient
    double detdF = dF.determinant();
	Matrix3 pF = mptr->GetDeformationGradientMatrix();
    Matrix3 F = dR*pF;
    if(np==THREED_MPM)
        F.Scale(pow(detdF,1./3.));
    else
        F.Scale2D(sqrt(detdF));
	
	// new deformation matrix with volume change onle
    mptr->SetDeformationGradientMatrix(F);
#else

	// Update total deformation gradient, and calculate trial B
	double detdF = IncrementDeformation(mptr,du,NULL,np);
#endif
    
    // Get new J and save result on the particle
	double J = detdF * mptr->GetHistoryDble(J_history);
    mptr->SetHistoryDble(J_history,J);
    
    // account for residual stresses
	double dresStretch,resStretch = GetResidualStretch(mptr,dresStretch,res);
	double Jres = resStretch*resStretch*resStretch;
    double Jeff = J/Jres;

    // new Kirchhoff pressure (over rho0) from Tait equation
	double p0=mptr->GetPressure();
    double pressure = J*TAIT_C*Ksp*(exp((1.-Jeff)/TAIT_C)-1.);
    mptr->SetPressure(pressure);
    
	// incremental energy per unit mass - dilational part
    double avgP = 0.5*(p0+pressure);
    double delV = 1. - 1./detdF;
    double workEnergy = -avgP*delV;
    
	// incremental residual energy per unit mass
	double delVres = 1. - 1./(dresStretch*dresStretch*dresStretch);
	double resEnergy = -avgP*delVres;
	
    // viscosity term = 2 eta (0.5(grad v) + 0.5*(grad V)^T - (1/3) tr(grad v) I)
    // (i.e., divatoric part of the symmetric strain tensor, 2 is for conversion to engineering shear strain)
    Matrix3 shear;
    double c[3][3];
    c[0][0] = (2.*du(0,0)-du(1,1)-du(2,2))/3.;
    c[1][1] = (2.*du(1,1)-du(0,0)-du(2,2))/3.;
    c[2][2] = (2.*du(2,2)-du(0,0)-du(1,1))/3.;
    c[0][1] = 0.5*(du(0,1)+du(1,0));
    c[1][0] = c[0][1];
    if(np==THREED_MPM)
    {   c[0][2] = 0.5*(du(0,2)+du(2,0));
        c[2][0] = c[0][2];
        c[1][2] = 0.5*(du(1,2)+du(2,1));
        c[2][1] = c[1][2];
        shear.set(c);
    }
    else
        shear.set(c[0][0],c[0][1],c[1][0],c[1][1],c[2][2]);
    
    // Get Kirchoff shear stress (over rho0)
    shear.Scale(J*TwoEtasp/delTime);
    
    // update deviatoric stress
	Tensor *sp=mptr->GetStressTensor();
    sp->xx = shear(0,0);
    sp->yy = shear(1,1);
    sp->zz = shear(2,2);
    sp->xy = shear(0,1);
    if(np==THREED_MPM)
    {   sp->xz = shear(0,2);
        sp->yz = shear(1,2);
    }
    
    // shear work per unit mass = tau.du = tau.tau*delTime/TwoEtasp
    double shearWork = sp->xx*sp->xx + sp->yy*sp->yy + sp->zz*sp->zz + 2.*sp->xy*sp->xy;
    if(np==THREED_MPM) shearWork += 2.*(sp->xz*sp->xz + sp->yz*sp->yz);
    shearWork *= delTime/TwoEtasp;
    mptr->AddWorkEnergyAndResidualEnergy(workEnergy+shearWork,resEnergy);
    
    // particle isentropic temperature increment dT/T = - J (K/K0) gamma0 Delta(V)/V
    // Delta(V)/V = 1. - 1/detdF (total volume)
	double Kratio = Jeff*(1.+pressure/(TAIT_C*Ksp));
	double dTq0 = -J*Kratio*gamma0*mptr->pPreviousTemperature*delV;
    
    // heat energy is Cv (dT - dTq0) -dPhi
	// Here do Cv (dT - dTq0)
    // dPhi = shearWork is lost due to shear term
    IncrementHeatEnergy(mptr,res->dT,dTq0,shearWork);
}

// Entry point for large rotation
void Elastic::LRConstitutiveLaw(MPMBase *mptr,Matrix3 du,double delTime,int np,void *properties,ResidualStrains *res) const
{
	// get previous deformation gradient
	Matrix3 pFnm1 = mptr->GetDeformationGradientMatrix();
	
    // get incremental deformation gradient and decompose it
	const Matrix3 dF = du.Exponential(incrementalDefGradTerms);
    Matrix3 dR;
    Matrix3 dU = dF.RightDecompose(&dR,NULL);

	// get pinitial rotation R0
	Matrix3 R0 = mptr->GetInitialRotation();
	
	// get previous rotation and stretch
    Matrix3 Rnm1;
	Matrix3 Unm1 = pFnm1.RightDecompose(&Rnm1,NULL);
	
	// get strain increments de = R0T.[(Rnm1T.dU.Rnm1 - I).Unm1].R0
	Matrix3 dUrot = dU.RTMR(Rnm1);
	dUrot(0,0) -= 1.;
	dUrot(1,1) -= 1.;
	dUrot(2,2) -= 1.;
	dUrot *= Unm1;
	
	// apply initial rotation to get strain increment in the material coordinates
	Matrix3 de = dUrot.RTMR(R0);
	Matrix3 Rtotnm1M3 = Rnm1*R0;
	Matrix3 *Rtotnm1 = &Rtotnm1M3;
	
	// get total rotation
	Matrix3 Rtot = dR*Rtotnm1M3;
	if(np==THREED_MPM) mptr->SetRtot(Rtot);
	
	// Update total deformation gradient
	Matrix3 pF = dF*pFnm1;
	mptr->SetDeformationGradientMatrix(pF);
	
	// cast pointer to material-specific data
	ElasticProperties *p = GetElasticPropertiesPointer(properties);

    // residual strains (thermal and moisture) in material axes
	double exxr,eyyr,ezzr;
	if(np==THREED_MPM)
	{	exxr = p->alpha[0]*res->dT;
		eyyr = p->alpha[1]*res->dT;
		ezzr = p->alpha[2]*res->dT;
		if(DiffusionTask::active)
		{	exxr += p->beta[0]*res->dC;
			eyyr += p->beta[1]*res->dC;
			ezzr += p->beta[2]*res->dC;
		}
	}
	else
	{	exxr = p->alpha[1]*res->dT;
		eyyr = p->alpha[2]*res->dT;
		ezzr = p->alpha[4]*res->dT;
		if(DiffusionTask::active)
		{	exxr += p->beta[1]*res->dC;
			eyyr += p->beta[2]*res->dC;
			ezzr += p->beta[4]*res->dC;
		}
	}
	Matrix3 er = Matrix3(exxr,0.,0.,eyyr,ezzr);
	
	// finish up
	LRElasticConstitutiveLaw(mptr,de,er,Rtot,dR,Rtotnm1,np,properties,res);
}