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

void AxisViewportRect(ViewExp *vpt, const Matrix3 &tm, float length, Rect *rect)
	{
	Matrix3 tmn = tm;
	float zoom;
	IPoint3 wpt;
	Point3 pt;
	GraphicsWindow *gw = vpt->getGW();

	// Get width of viewport in world units:  --DS
	zoom = vpt->GetScreenScaleFactor(tmn.GetTrans())*ZFACT;
	
	tmn.Scale( Point3(zoom,zoom,zoom) );
	gw->setTransform( tmn );	
	pt = Point3(0.0f, 0.0f, 0.0f);
	gw->wTransPoint( &pt, &wpt );
	rect->left = rect->right  = wpt.x;
	rect->top  = rect->bottom = wpt.y;

	AxisRect( gw, Point3(length,0.0f,0.0f),rect );	
	AxisRect( gw, Point3(0.0f,length,0.0f),rect );	
	AxisRect( gw, Point3(0.0f,0.0f,length),rect );	

	rect->right  += 2;
	rect->bottom += 2;
	rect->left   -= 2;
	rect->top    -= 2;
	}

void ProtHelpObject::GetMat(TimeValue t, INode* inode, ViewExp& vpt, Matrix3& tm) 
{
   if ( ! vpt.IsAlive() )
	{
		tm.Zero();
		return;
	}
	 tm = inode->GetObjectTM(t);
   tm.NoScale();
   float scaleFactor = vpt.NonScalingObjectSize() * vpt.GetVPWorldWidth(tm.GetTrans()) / 360.0f;
   tm.Scale(Point3(scaleFactor,scaleFactor,scaleFactor));
}

void VRayCamera::getTM(TimeValue t, INode *node, ViewExp *vpt, Matrix3 &tm) {
  tm=node->GetObjectTM(t);

  AffineParts ap;
  decomp_affine(tm, &ap);
  tm.IdentityMatrix();
  tm.SetRotate(ap.q);
  tm.SetTrans(ap.t);

  float scaleFactor=vpt->NonScalingObjectSize()*vpt->GetVPWorldWidth(tm.GetTrans())/360.0f;
  tm.Scale(Point3(scaleFactor,scaleFactor,scaleFactor));
}

void SymmetryMod::MirrorTriObject (Mesh & mesh, int axis, Matrix3 & tm, Matrix3 & itm,  int normalMapChannel) {
	// Create scaling matrix for mirroring on selected axis:
	Point3 scale(1,1,1);
	scale[axis] = -1.0f;
	itm.Scale(scale,TRUE);

	// Hang on to a copy of the incoming face selection:
	BitArray inputFaceSel = mesh.faceSel;

	// Make the mirror copy of the entire mesh:
	int oldnumv = mesh.numVerts;
	int oldnumf = mesh.numFaces;

	int oldNumNormals = 0;
	if (normalMapChannel != INVALID_NORMALMAPCHANNEL)	
	{
		MeshMap& map = mesh.Map(normalMapChannel);
		oldNumNormals = map.vnum;
	}

	BitArray fset;
	fset.SetSize (oldnumf);
	fset.SetAll ();
	mesh.CloneFaces (fset);	// Clears selection on originals, sets it on new faces.

	// Transform the cloned vertices to their mirror images:
	for (int i=oldnumv; i<mesh.numVerts; i++) {
		mesh.verts[i] = (mesh.verts[i]*itm)*tm;
	}

	// Restore selection of input faces:
	for (int i=0; i<oldnumf; i++) mesh.faceSel.Set (i, inputFaceSel[i]);
	// Flip over new faces and select to match input:
	for (int i=oldnumf; i<mesh.numFaces; i++) {
		mesh.FlipNormal (i);
		mesh.faceSel.Set (i, inputFaceSel[i-oldnumf]);
	}

	//flip and specified normals/faces
	if (normalMapChannel != INVALID_NORMALMAPCHANNEL)	
	{
		MeshMap& map = mesh.Map(normalMapChannel);
		int numNormals = map.vnum;
		
		Matrix3 mirrorTM = itm*tm;		
		for (int i = oldNumNormals; i < numNormals; i++)
		{
			Point3 n = map.tv[i];
			n = VectorTransform(n,mirrorTM);
			map.tv[i] = n;
		}
	}
}

// Update Strains for this particle
// Velocities for all fields are present on the nodes
// matRef is the material and properties have been loaded, matFld is the material field
void MatPointAS::UpdateStrain(double strainTime,int secondPass,int np,void *props,int matFld)
{
	int i,numnds,nds[maxShapeNodes];
    double fn[maxShapeNodes],xDeriv[maxShapeNodes],yDeriv[maxShapeNodes],zDeriv[maxShapeNodes];
	Vector vel;
    Matrix3 dv;
	
	// don't need to zero zDeriv because always set in axisymmetric shape functions
    
	// find shape functions and derviatives
	const ElementBase *elemRef = theElements[ElemID()];
	elemRef->GetShapeGradients(&numnds,fn,nds,xDeriv,yDeriv,zDeriv,this);
    
    // Find strain rates at particle from current grid velocities
	//   and using the velocity field for that particle and each node and the right material
    // In axisymmetric x->r, y->z, and z->hoop
    for(i=1;i<=numnds;i++)
	{	vel=nd[nds[i]]->GetVelocity((short)vfld[i],matFld);
        dv += Matrix3(vel.x*xDeriv[i],vel.x*yDeriv[i],vel.y*xDeriv[i],vel.y*yDeriv[i],vel.x*zDeriv[i]);
    }
    
    // save velocity gradient (if needed for J integral calculation)
    SetVelocityGradient(dv(0,0),dv(1,0),dv(0,1),dv(1,0),secondPass);
    
    // convert to strain increments
    // e.g., now dvrr = dvr/dr * dt = d/dr(du/dt) * dt = d/dt(du/dr) * dt = du/dr)
    dv.Scale(strainTime);
    
	// find effective particle transport properties from grid results
	ResidualStrains res;
	res.dT = 0;
	res.dC = 0.;
	if(!ConductionTask::active)
	{	res.dT = pTemperature-pPreviousTemperature;
		pPreviousTemperature = pTemperature;
	}
	else
	{	for(i=1;i<=numnds;i++)
		res.dT += conduction->IncrementValueExtrap(nd[nds[i]],fn[i]);
		res.dT = conduction->GetDeltaValue(this,res.dT);
	}
	if(DiffusionTask::active)
	{	for(i=1;i<=numnds;i++)
		res.dC += diffusion->IncrementValueExtrap(nd[nds[i]],fn[i]);
		res.dC = diffusion->GetDeltaValue(this,res.dC);
	}
	
	// pass on to material class to handle
	PerformConstitutiveLaw(dv,strainTime,np,props,&res);
}

void SymmetryMod::MirrorPolyObject (MNMesh & mesh, int axis, Matrix3 & tm, Matrix3 & itm) {
	// Create scaling matrix for mirroring on selected axis:
	Point3 scale(1,1,1);
	scale[axis] = -1.0f;
	itm.Scale(scale,TRUE);

	// Make the mirror copy of the entire mesh:
	int oldnumv = mesh.numv;
	for (int i=0; i<mesh.numf; i++) mesh.f[i].SetFlag (MN_USER);
	mesh.CloneFaces (MN_USER, true);

	// Transform the vertices to their mirror images:
	for (int i=oldnumv; i<mesh.numv; i++) mesh.v[i].p = (itm*mesh.v[i].p)*tm;

	// Flip over faces, edges:
	mesh.FlipElementNormals (MN_USER);	// flag should now be set only on clones.

	DbgAssert (mesh.CheckAllData ());
}

void DrawAxis(ViewExp *vpt, const Matrix3 &tm, float length, BOOL screenSize)
	{
	Matrix3 tmn = tm;
	float zoom;

	// Get width of viewport in world units:  --DS
	zoom = vpt->GetScreenScaleFactor(tmn.GetTrans())*ZFACT;
	
	if (screenSize) {
		tmn.Scale( Point3(zoom,zoom,zoom) );
		}
	vpt->getGW()->setTransform( tmn );		

	Text( vpt, _T("x"), Point3(length,0.0f,0.0f) ); 
	DrawAnAxis( vpt, Point3(length,0.0f,0.0f) );	
	
	Text( vpt, _T("y"), Point3(0.0f,length,0.0f) ); 
	DrawAnAxis( vpt, Point3(0.0f,length,0.0f) );	
	
	Text( vpt, _T("z"), Point3(0.0f,0.0f,length) ); 
	DrawAnAxis( vpt, Point3(0.0f,0.0f,length) );
	}

bool Exporter::splitMesh(INode *node, Mesh& mesh, FaceGroups &grps, TimeValue t, vector<Color4>& vertColors, bool noSplit)
{
	Mtl* nodeMtl = node->GetMtl();
	Matrix3 tm = node->GetObjTMAfterWSM(t);

	// Order of the vertices. Get 'em counter clockwise if the objects is
	// negatively scaled.
	int vi[3];
	if (TMNegParity(tm)) {
		vi[0] = 2; vi[1] = 1; vi[2] = 0;
	} else {
		vi[0] = 0; vi[1] = 1; vi[2] = 2;
	}

	Matrix3 flip;
	flip.IdentityMatrix();
	flip.Scale(Point3(1, -1, 1));

   int nv = mesh.getNumVerts();
   int nf = mesh.getNumFaces();

	if (noSplit)
	{
		int nv = mesh.getNumVerts();
		int nf = mesh.getNumFaces();
		// Dont split the mesh at all.  For debugging purposes.
		FaceGroup& grp = grps[0];
		grp.vidx.resize(nv, -1);
		grp.verts.resize(nv);
		grp.faces.resize(nf);
		grp.uvs.resize(nv);
		grp.vnorms.resize(nv);
      grp.fidx.resize(nf);

		Matrix3 texm;
		getTextureMatrix(texm, getMaterial(node, 0));
		texm *= flip;

		for (int face=0; face<nf; ++face) {
         grp.fidx[face] = face;
			for (int vi=0; vi<3; ++vi) {
				int idx = mesh.faces[face].getVert(vi);
				grp.faces[face][vi] = idx;

				// Calculate normal
				Point3 norm;
#if VERSION_3DSMAX <= ((5000<<16)+(15<<8)+0) // Version 5
				norm = getVertexNormal(&mesh, face, mesh.getRVertPtr(idx));
#else
				MeshNormalSpec *specNorms = mesh.GetSpecifiedNormals ();
				if (NULL != specNorms && specNorms->GetNumNormals() != 0)
					norm = specNorms->GetNormal(face, vi);
				else
					norm = getVertexNormal(&mesh, face, mesh.getRVertPtr(idx));
#endif
				Point3 uv;
				if (mesh.tVerts && mesh.tvFace) {
					uv = mesh.tVerts[ mesh.tvFace[ face ].t[ vi ]] * texm;
					uv.y += 1.0f;
				}

				if (grp.vidx[idx] == idx){
					ASSERT(grp.verts[idx] == TOVECTOR3(mesh.getVert(idx)));
					//ASSERT(vg.norm == norm);
					//Point3 uv = mesh.getTVert(idx);
					//if (mesh.getNumTVerts() > 0)
					//{
					//	ASSERT(grp.uvs[idx].u == uv.x && grp.uvs[idx].v == uv.y);
					//}
				} else {
					grp.vidx[idx] = idx;
					grp.verts[idx] = TOVECTOR3(mesh.getVert(idx));
					//grp.uvs[idx].u = uv.x;
					//grp.uvs[idx].v = uv.y;
					grp.vnorms[idx] = TOVECTOR3(norm);
				}
			}
		}
		for (int i=0; i<nv; ++i) {
			ASSERT(grp.vidx[i] != -1);
		}
	}
	else
	{
		int face, numSubMtls = nodeMtl?nodeMtl->NumSubMtls():0;
		for (face=0; face<mesh.getNumFaces(); face++) 
		{
			int mtlID = (numSubMtls!=0) ? (mesh.faces[face].getMatID() % numSubMtls) : 0;
         Mtl *mtl = getMaterial(node, mtlID);
			Matrix3 texm;
			getTextureMatrix(texm, mtl);
			texm *= flip;

         FaceGroup& grp = grps[mtlID];

		 if (grp.uvMapping.size() == 0) // Only needs to be done once per face group
		 {
			 int nmaps = 0;
			 int nmapsStart = max(1, mesh.getNumMaps() - (mesh.mapSupport(0) ? 1 : 0)); // Omit vertex color map.
			 for (int ii = 1; ii <= nmapsStart; ii++) // Winnow out the unsupported maps.
//.........这里部分代码省略.........

// 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];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	// make verts
	for(ix=0; ix<segs; ix++) {
		sinang = (float)sin(ang);
		cosang = (float)cos(ang);
		ang2 = rotation + twist * float(ix+1)/float(segs);
		for (jx = 0; jx<sides; jx++) {
			sinang2 = (float)sin(ang2);
			cosang2 = (float)cos(ang2);
			rt = radius+radius2*cosang2;
			p.x = rt*cosang;
			p.y = -rt*sinang;
			p.z = radius2*sinang2;	
			mesh.setVert(nv++, p);
			ang2 += delta2;
			}	
		ang += delta;
		}
	
	if (doPie) {
		p.x = radius * (float)cos(startAng);
		p.y = -radius * (float)sin(startAng);
		p.z = 0.0f;
		mesh.setVert(nv++, p);

		ang -= delta;
		p.x = radius * (float)cos(ang);
		p.y = -radius * (float)sin(ang);
		p.z = 0.0f;
		mesh.setVert(nv++, p);
		}
	
	// Make faces

    BOOL usePhysUVs = GetUsePhysicalScaleUVs();
    BitArray startSliceFaces;
    BitArray endSliceFaces;

    if (usePhysUVs) {
        startSliceFaces.SetSize(mesh.numFaces);
        endSliceFaces.SetSize(mesh.numFaces);
    }

	/* Make midsection */
	for(ix=0; ix<segs-doPie; ++ix) {
		jx=ix*sides;
		for (kx=0; kx<sides; ++kx) {
			na = jx+kx;
			nb = (ix==(segs-1))?kx:na+sides;
			nd = (kx==(sides-1))? jx : na+1;
			nc = nb+nd-na;

			DWORD grp = 0;
			if (smooth==SMOOTH_SIDES) {
				if (kx==sides-1 && (sides&1)) {
					grp = (1<<2);
				} else {
					grp = (kx&1) ? (1<<0) : (1<<1);
					}
			} else 
			if (smooth==SMOOTH_STRIPES) {
				if (ix==segs-1 && (segs&1)) {
					grp = (1<<2);
				} else {
					grp = (ix&1) ? (1<<0) : (1<<1);
					}
			} else 

Object*
BuildNURBSPyramid(float width, float depth, float height, int genUVs)
{
	int pyramid_faces[5][4] = { {0, 1, 2, 3}, // bottom
							{2, 3, 4, 4}, // back
							{1, 0, 4, 4}, // front
							{3, 1, 4, 4}, // left
							{0, 2, 4, 4}};// right
	Point3 pyramid_verts[5] = { Point3(-0.5, -0.5, 0.0),
							Point3( 0.5, -0.5, 0.0),
							Point3(-0.5,  0.5, 0.0),
							Point3( 0.5,  0.5, 0.0),
							Point3( 0.0,  0.0, 1.0)};

	NURBSSet nset;

	for (int face = 0; face < 5; face++) {
		Point3 bl = pyramid_verts[pyramid_faces[face][0]];
		Point3 br = pyramid_verts[pyramid_faces[face][1]];
		Point3 tl = pyramid_verts[pyramid_faces[face][2]];
		Point3 tr = pyramid_verts[pyramid_faces[face][3]];

		Matrix3 size;
		size.IdentityMatrix();
		Point3 lwh(width, depth, height);
		size.Scale(lwh);

		bl = bl * size;
		br = br * size;
		tl = tl * size;
		tr = tr * size;

		NURBSCVSurface *surf = new NURBSCVSurface();
		nset.AppendObject(surf);
		surf->SetUOrder(4);
		surf->SetVOrder(4);
		surf->SetNumCVs(4, 4);
		surf->SetNumUKnots(8);
		surf->SetNumVKnots(8);

		Point3 top, bot;
		for (int r = 0; r < 4; r++) {
			top = tl + (((float)r/3.0f) * (tr - tl));
			bot = bl + (((float)r/3.0f) * (br - bl));
			for (int c = 0; c < 4; c++) {
				NURBSControlVertex ncv;
				ncv.SetPosition(0, bot + (((float)c/3.0f) * (top - bot)));
				ncv.SetWeight(0, 1.0f);
				surf->SetCV(r, c, ncv);
			}
		}

		for (int k = 0; k < 4; k++) {
			surf->SetUKnot(k, 0.0);
			surf->SetVKnot(k, 0.0);
			surf->SetUKnot(k + 4, 1.0);
			surf->SetVKnot(k + 4, 1.0);
		}

		surf->Renderable(TRUE);
		surf->SetGenerateUVs(genUVs);
		if (height > 0.0f)
			surf->FlipNormals(TRUE);
		else
			surf->FlipNormals(FALSE);

		switch(face) {
		case 0: // bottom
			surf->SetTextureUVs(0, 0, Point2(1.0f, 0.0f));
			surf->SetTextureUVs(0, 1, Point2(0.0f, 0.0f));
			surf->SetTextureUVs(0, 2, Point2(1.0f, 1.0f));
			surf->SetTextureUVs(0, 3, Point2(0.0f, 1.0f));
			break;
		default: // sides
			surf->SetTextureUVs(0, 0, Point2(0.5f, 1.0f));
			surf->SetTextureUVs(0, 1, Point2(0.5f, 1.0f));
			surf->SetTextureUVs(0, 2, Point2(0.0f, 0.0f));
			surf->SetTextureUVs(0, 3, Point2(1.0f, 0.0f));
			break;
		}

		TCHAR bname[80];
		_stprintf(bname, _T("%s%02d"), GetString(IDS_CT_SURF), face);
		surf->SetName(bname);
	}

#define F(s1, s2, s1r, s1c, s2r, s2c) \
	fuse.mSurf1 = (s1); \
	fuse.mSurf2 = (s2); \
	fuse.mRow1 = (s1r); \
	fuse.mCol1 = (s1c); \
	fuse.mRow2 = (s2r); \
	fuse.mCol2 = (s2c); \
	nset.mSurfFuse.Append(1, &fuse);

	NURBSFuseSurfaceCV fuse;
	// Fuse the degenerate peaks
	for (int i = 1; i < 5; i++) {
		for (int j = 1; j < 4; j++) {
			F(i, i, 0, 3, j, 3);
//.........这里部分代码省略.........

// 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);
}