Method updateActiveSurface

SRC/material/nD/soil/MultiYieldSurfaceClay.cpp:1756–1827

void MultiYieldSurfaceClay::updateActiveSurface(void)
{
  int numOfSurfaces = numOfSurfacesx[matN];

  if (activeSurfaceNum == numOfSurfaces) return;

	double A, B, C, X;
	static T2Vector direction;
	static Vector t1(6);
	static Vector t2(6);
//	static Vector temp(6);
	static Vector center(6);
	center = theSurfaces[activeSurfaceNum].center();
	double size = theSurfaces[activeSurfaceNum].size();
	static Vector outcenter(6);
	outcenter= theSurfaces[activeSurfaceNum+1].center();
	double outsize = theSurfaces[activeSurfaceNum+1].size();


	//t1 = trialStress.deviator() - center;
	//t2 = center - outcenter;
	t1 = trialStress.deviator();
	t1 -= center;
	t2 = center;
	t2 -= outcenter;

	A = t1 && t1;
	B = 2. * (t1 && t2);
	C = (t2 && t2) - 2./3.* outsize * outsize;
	X = secondOrderEqn(A,B,C,0);
	if ( fabs(X-1.) < LOW_LIMIT ) X = 1.;
	if (X < 1.){
	  opserr << "FATAL:MultiYieldSurfaceClay::updateActiveSurface(): error in Direction of surface motion." 
	       << endln; 
	  exit(-1);
	}

	//temp = (t1 * X + center) * (1. - size / outsize) - (center - outcenter * size / outsize);
	temp = center;
	temp.addVector(1.0, t1, X);
	temp *= (1.0 - size/outsize);
	t2 = center;
	t2.addVector(1.0, outcenter, -size/outsize);
	temp -= t2;

	direction.setData(temp);

	if (direction.deviatorLength() < LOW_LIMIT) return;

	temp = direction.deviator();  
	A = temp && temp;
	B = - 2 * (t1 && temp);
	if (fabs(B) < LOW_LIMIT) B = 0.; 
	C = (t1 && t1) - 2./3.* size * size;
	if ( fabs(C) < LOW_LIMIT || fabs(C)/(t1 && t1) < LOW_LIMIT ) return;

// guquan 2005 Apr. 21
//	opserr<<"Low_limit is " << LOW_LIMIT << endln;