C++ StructuralMaterialStatus::giveTempStrainVector方法代码示例

void
HyperElasticMaterial :: give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
{
    double J2, c11, c22, c33, c12, c13, c23, A, B;
    FloatMatrix C(3, 3);
    FloatMatrix invC;

    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );

    C.at(1, 1) = 1. + 2. * status->giveTempStrainVector().at(1);
    C.at(2, 2) = 1. + 2. * status->giveTempStrainVector().at(2);
    C.at(3, 3) = 1. + 2. * status->giveTempStrainVector().at(3);
    C.at(2, 3) = C.at(3, 2) = status->giveTempStrainVector().at(4);
    C.at(1, 3) = C.at(3, 1) = status->giveTempStrainVector().at(5);
    C.at(1, 2) = C.at(2, 1) = status->giveTempStrainVector().at(6);

    invC.beInverseOf(C);
    J2 = C.giveDeterminant();

    c11 = invC.at(1, 1);
    c22 = invC.at(2, 2);
    c33 = invC.at(3, 3);
    c12 = invC.at(1, 2);
    c13 = invC.at(1, 3);
    c23 = invC.at(2, 3);

    A = ( K - 2. / 3. * G ) * J2;
    B = -( K - 2. / 3. * G ) * ( J2 - 1. ) + 2. * G;

    answer.resize(6, 6);

    answer.at(1, 1) = ( A + B ) * c11 * c11;
    answer.at(2, 2) = ( A + B ) * c22 * c22;
    answer.at(3, 3) = ( A + B ) * c33 * c33;
    answer.at(4, 4) = A * c23 * c23 + B / 2. * ( c22 * c33 + c23 * c23 );
    answer.at(5, 5) = A * c13 * c13 + B / 2. * ( c11 * c33 + c13 * c13 );
    answer.at(6, 6) = A * c12 * c12 + B / 2. * ( c11 * c22 + c12 * c12 );
    answer.at(1, 2) = answer.at(2, 1) = A * c11 * c22 + B * c12 * c12;
    answer.at(1, 3) = answer.at(3, 1) = A * c11 * c33 + B * c13 * c13;
    answer.at(1, 4) = answer.at(4, 1) = A * c11 * c23 + B * c12 * c13;
    answer.at(1, 5) = answer.at(5, 1) = A * c11 * c13 + B * c11 * c13;
    answer.at(1, 6) = answer.at(6, 1) = A * c11 * c12 + B * c11 * c12;
    answer.at(2, 3) = answer.at(3, 2) = A * c22 * c33 + B * c23 * c23;
    answer.at(2, 4) = answer.at(4, 2) = A * c22 * c23 + B * c22 * c23;
    answer.at(2, 5) = answer.at(5, 2) = A * c22 * c13 + B * c12 * c23;
    answer.at(2, 6) = answer.at(6, 2) = A * c22 * c12 + B * c22 * c12;
    answer.at(3, 4) = answer.at(4, 3) = A * c33 * c23 + B * c33 * c23;
    answer.at(3, 5) = answer.at(5, 3) = A * c33 * c13 + B * c33 * c13;
    answer.at(3, 6) = answer.at(6, 3) = A * c33 * c12 + B * c13 * c23;
    answer.at(4, 5) = answer.at(5, 4) = A * c23 * c13 + B / 2. * ( c12 * c33 + c13 * c23 );
    answer.at(4, 6) = answer.at(6, 4) = A * c23 * c12 + B / 2. * ( c12 * c23 + c22 * c13 );
    answer.at(5, 6) = answer.at(6, 5) = A * c13 * c12 + B / 2. * ( c11 * c23 + c12 * c13 );
}

double
PhaseFieldElement :: computeFreeEnergy(GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );
    FloatArray strain, stress;
    stress = matStat->giveTempStressVector();
    strain = matStat->giveTempStrainVector();
    return 0.5 * stress.dotProduct( strain );
}

std::vector<std::unique_ptr<EnrichmentItem>> NCPrincipalStrain::nucleateEnrichmentItems() {


	SpatialLocalizer *octree = this->mpDomain->giveSpatialLocalizer();
	XfemManager *xMan = mpDomain->giveXfemManager();

	std::vector<std::unique_ptr<EnrichmentItem>> eiList;

	// Center coordinates of newly inserted cracks
	std::vector<FloatArray> center_coord_inserted_cracks;

	// Loop over all elements and all bulk GP.
	for(auto &el : mpDomain->giveElements() ) {

		int numIR = el->giveNumberOfIntegrationRules();

		int csNum = el->giveCrossSection()->giveNumber();

		if(csNum == mCrossSectionInd) {

			for(int irInd = 0; irInd < numIR; irInd++) {
				IntegrationRule *ir = el->giveIntegrationRule(irInd);


				int numGP = ir->giveNumberOfIntegrationPoints();

				for(int gpInd = 0; gpInd < numGP; gpInd++) {
					GaussPoint *gp = ir->getIntegrationPoint(gpInd);


						StructuralMaterialStatus *ms = dynamic_cast<StructuralMaterialStatus*>(gp->giveMaterialStatus());

						if(ms != NULL) {

							const FloatArray &strain = ms->giveTempStrainVector();

							FloatArray principalVals;
							FloatMatrix principalDirs;
							StructuralMaterial::computePrincipalValDir(principalVals, principalDirs, strain, principal_strain);

							if(principalVals[0] > mStrainThreshold) {

								FloatArray crackNormal;
								crackNormal.beColumnOf(principalDirs, 1);
		//						printf("crackNormal: "); crackNormal.printYourself();

								FloatArray crackTangent = {-crackNormal(1), crackNormal(0)};
								crackTangent.normalize();
		//						printf("crackTangent: "); crackTangent.printYourself();


								// Create geometry
								FloatArray pc = {gp->giveGlobalCoordinates()(0), gp->giveGlobalCoordinates()(1)};
		//						printf("Global coord: "); pc.printYourself();


								FloatArray ps = pc;
								ps.add(-0.5*mInitialCrackLength, crackTangent);

								FloatArray pe = pc;
								pe.add(0.5*mInitialCrackLength, crackTangent);

								if(mCutOneEl) {
									// If desired, ensure that the crack cuts exactly one element.
									Line line(ps, pe);
									std::vector<FloatArray> intersecPoints;
		//							line.computeIntersectionPoints(el.get(), intersecPoints);


									if(intersecPoints.size() == 2) {
										ps = std::move(intersecPoints[0]);
										pe = std::move(intersecPoints[1]);
									}
									else {
										OOFEM_ERROR("intersecPoints.size() != 2")
									}
								}

								FloatArray points = {ps(0), ps(1), pc(0), pc(1), pe(0), pe(1)};


								// Check if nucleation is allowed, by checking for already existing cracks close to the GP.
								// Idea: Nucleation is not allowed if we are within an enriched element. In this way, branching is not
								// completely prohibited, but we avoid initiating multiple similar cracks.
								bool insertionAllowed = true;

								Element *el_s = octree->giveElementContainingPoint(ps);
								if(el_s) {
									if( xMan->isElementEnriched(el_s) ) {
										insertionAllowed = false;
									}
								}

								Element *el_c = octree->giveElementContainingPoint(pc);
								if(el_c) {
									if( xMan->isElementEnriched(el_c) ) {
										insertionAllowed = false;
									}
								}

//.........这里部分代码省略.........