Method compute

SRC/material/nD/ASDConcrete3DMaterial.cpp:781–852

}

int ASDConcrete3DMaterial::StressDecomposition::compute(const Vector& S, double cdf)
{
	// copy S in V
	V(0, 0) = S(0);
	V(1, 1) = S(1);
	V(2, 2) = S(2);
	V(0, 1) = V(1, 0) = S(3);
	V(1, 2) = V(2, 1) = S(4);
	V(0, 2) = V(2, 0) = S(5);

	// solve
	if (Eigen3(Si, V) < 0)
		return EC_Eigen_Error;

	// construct matrices PT and PC
	static Matrix pjj(6, 6);
	PT.Zero();
	PC.Zero();

	// method 1: compute PT = sum(H(sj)*pjj) --- PC = I - PT
	//for (int j = 0; j < 3; j++) {
	//	double hsj = Heavyside(Si(j));
	//	if (hsj > 0.0) {
	//		computePjj(V, j, pjj);
	//		PT.addMatrix(1.0, pjj, hsj);
	//	}
	//}
	//for (int i = 0; i < 6; ++i)
	//	PC(i, i) = 1.0;
	//PC.addMatrix(1.0, PT, -1.0);

	// method 2: compute PT = sum(H(sj)*pjj) --- PC = sum(H(-sj)*pjj)
	// then PO = I - PT - PC
	// finally PT += PO/2 --- PC += PO/2
	for (int j = 0; j < 3; j++) {
		computePjj(V, j, pjj);
		double hsj = Heavyside(Si(j));
		if(hsj > 0.0)
			PT.addMatrix(1.0, pjj, hsj);
		hsj = Heavyside(-Si(j));
		if (hsj > 0.0)
			PC.addMatrix(1.0, pjj, hsj);
	}

	static Matrix PO(6, 6); // PO = I - PT - PC
	PO.addMatrix(0.0, PT, -1.0);
	PO.addMatrix(1.0, PC, -1.0);
	for (int i = 0; i < 6; ++i)
		PO(i, i) += 1.0;
	PT.addMatrix(1.0, PO, 0.5);
	PC.addMatrix(1.0, PO, 0.5);

	// compute positive and negative parts of the stress
	ST.addMatrixVector(0.0, PT, S, 1.0);
	SC.addMatrixVector(0.0, PC, S, 1.0);

	// R factor
	R = 0.0;