WebSVN - OpenSees - Blame - Rev 458 - /trunk/SRC/material/nD/Template3Dep/Template3Dep.cpp

Rev	Author	Line No.	Line
12	jeremic	1	/*
		2	################################################################################
		3	# COPYRIGHT (C): :-)) #
		4	# PROJECT: Object Oriented Finite Element Program #
		5	# PURPOSE: General platform for elaso-plastic constitutive model #
		6	# implementation #
		7	# CLASS: Template3Dep (the base class for all material point) #
		8	# #
		9	# VERSION: #
		10	# LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) #
		11	# TARGET OS: DOS \|\| UNIX \|\| . . . #
		12	# DESIGNER(S): Boris Jeremic, Zhaohui Yang #
		13	# PROGRAMMER(S): Boris Jeremic, Zhaohui Yang #
		14	# #
		15	# #
		16	# DATE: 08-03-2000 #
		17	# UPDATE HISTORY: 09-12-2000 #
		18	# #
		19	# #
		20	# #
		21	# #
		22	# SHORT EXPLANATION: This file contains the class implementation for #
		23	# Template3Dep. #
		24	# #
		25	################################################################################
		26	*/
		27
		28	#ifndef Template3Dep_CPP
		29	#define Template3Dep_CPP
		30
13	jeremic	31	#define ITMAX 30
		32	#define MAX_STEP_COUNT 30
130	jeremic	33	#define NUM_OF_SUB_INCR 10
12	jeremic	34	//#include <string.h>
		35
		36	#include "Template3Dep.h"
		37	#include <Channel.h>
		38	#include <G3Globals.h>
		39
		40	#include <DP_YS.h>
		41	#include <DP_PS.h>
		42	#include <EPState.h>
		43
		44	//================================================================================
		45	// Constructor
		46	//================================================================================
		47
13	jeremic	48	Template3Dep::Template3Dep( int tag ,
12	jeremic	49	YieldSurface *YS_ ,
		50	PotentialSurface *PS_ ,
		51	EPState *EPS_ ,
		52	EvolutionLaw_S *ELS1_ ,
		53	EvolutionLaw_S *ELS2_ ,
		54	EvolutionLaw_S *ELS3_ ,
		55	EvolutionLaw_S *ELS4_ ,
		56	EvolutionLaw_T *ELT1_ ,
		57	EvolutionLaw_T *ELT2_ ,
		58	EvolutionLaw_T *ELT3_ ,
		59	EvolutionLaw_T *ELT4_ )
		60	:NDMaterial(tag, ND_TAG_Template3Dep)
		61	{
		62	if ( YS_ )
		63	YS = YS_->newObj();
		64	else
		65	{
		66	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		67	exit(-1);
		68	}
		69
		70	if ( PS_ )
		71	PS = PS_->newObj();
		72	else
		73	{
		74	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		75	exit(-1);
		76	}
		77
		78	if ( EPS_ )
		79	EPS = EPS_->newObj();
		80	else
		81	{
		82	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		83	exit(-1);
		84	}
		85
		86	// Evolution laws
		87	if ( ELS1_ )
		88	ELS1 = ELS1_->newObj();
		89	else
		90	ELS1 = 0;
		91
		92	if ( ELS2_ )
		93	ELS2 = ELS2_->newObj();
		94	else
		95	ELS2 = 0;
		96
		97	if ( ELS3_ )
		98	ELS3 = ELS3_->newObj();
		99	else
		100	ELS3 = 0;
		101
		102	if ( ELS4_ )
		103	ELS4 = ELS4_->newObj();
		104	else
		105	ELS4 = 0;
		106
		107	if ( ELT1_ )
		108	ELT1 = ELT1_->newObj();
		109	else
		110	ELT1 = 0;
		111
		112	if ( ELT2_ )
		113	ELT2 = ELT2_->newObj();
		114	else
		115	ELT2 = 0;
		116
		117	if ( ELT3_ )
		118	ELT3 = ELT3_->newObj();
		119	else
		120	ELT3 = 0;
		121
		122	if ( ELT4_ )
		123	ELT4 = ELT4_->newObj();
		124	else
		125	ELT4 = 0;
		126
130	jeremic	127	//Initialze Eep using E-elastic
		128	tensor E = ElasticStiffnessTensor();
		129	EPS->setEep( E );
		130
12	jeremic	131	}
		132
		133	//================================================================================
		134	// Constructor 0
		135	//================================================================================
13	jeremic	136	Template3Dep::Template3Dep( int tag ,
12	jeremic	137	YieldSurface *YS_ ,
		138	PotentialSurface *PS_ ,
		139	EPState *EPS_)
		140	:NDMaterial(tag, ND_TAG_Template3Dep)
		141	{
		142	if ( YS_ )
		143	YS = YS_->newObj();
		144	else
		145	{
		146	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		147	exit(-1);
		148	}
		149
		150	if ( PS_ )
		151	PS = PS_->newObj();
		152	else
		153	{
		154	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		155	exit(-1);
		156	}
		157
		158	if ( EPS_ )
		159	EPS = EPS_->newObj();
		160	else
		161	{
		162	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		163	exit(-1);
		164	}
		165
		166	// Evolution laws
		167	ELS1 = 0;
		168	ELS2 = 0;
		169	ELS3 = 0;
		170	ELS4 = 0;
		171	ELT1 = 0;
		172	ELT2 = 0;
		173	ELT3 = 0;
		174	ELT4 = 0;
		175	}
		176
		177
		178	//================================================================================
		179	// Constructor 1
		180	//================================================================================
		181
13	jeremic	182	Template3Dep::Template3Dep( int tag ,
12	jeremic	183	YieldSurface *YS_ ,
		184	PotentialSurface *PS_ ,
		185	EPState *EPS_,
		186	EvolutionLaw_S *ELS1_ )
		187	:NDMaterial(tag, ND_TAG_Template3Dep)
		188	{
		189
		190	if ( YS_ )
		191	YS = YS_->newObj();
		192	else
		193	{
		194	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		195	exit(-1);
		196	}
		197
		198	if ( PS_ )
		199	PS = PS_->newObj();
		200	else
		201	{
		202	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		203	exit(-1);
		204	}
		205
		206	if ( EPS_ )
		207	EPS = EPS_->newObj();
		208	else
		209	{
		210	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		211	exit(-1);
		212	}
		213
		214	// Evolution laws
		215	if ( ELS1_ )
		216	ELS1 = ELS1_->newObj();
		217	else
		218	ELS1 = 0;
		219
		220	ELS2 = 0;
		221	ELS3 = 0;
		222	ELS4 = 0;
		223	ELT1 = 0;
		224	ELT2 = 0;
		225	ELT3 = 0;
		226	ELT4 = 0;
		227	}
		228
		229
		230	//================================================================================
		231	// Constructor 2
		232	//================================================================================
		233
13	jeremic	234	Template3Dep::Template3Dep( int tag ,
12	jeremic	235	YieldSurface *YS_ ,
		236	PotentialSurface *PS_ ,
		237	EPState *EPS_,
		238	EvolutionLaw_T *ELT1_ )
		239	:NDMaterial(tag, ND_TAG_Template3Dep)
		240	{
		241
		242	if ( YS_ )
		243	YS = YS_->newObj();
		244	else
		245	{
		246	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		247	exit(-1);
		248	}
		249
		250	if ( PS_ )
		251	PS = PS_->newObj();
		252	else
		253	{
		254	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		255	exit(-1);
		256	}
		257
		258	if ( EPS_ )
		259	EPS = EPS_->newObj();
		260	else
		261	{
		262	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		263	exit(-1);
		264	}
		265
		266	// Evolution laws
		267	ELS1 = 0;
		268	ELS2 = 0;
		269	ELS3 = 0;
		270	ELS4 = 0;
		271
		272	if ( ELT1_ )
		273	ELT1 = ELT1_->newObj();
		274	else
		275	ELT1 = 0;
		276
		277	ELT2 = 0;
		278	ELT3 = 0;
		279	ELT4 = 0;
		280	}
		281
		282	//================================================================================
		283	// Constructor 3
		284	//================================================================================
		285
13	jeremic	286	Template3Dep::Template3Dep( int tag ,
12	jeremic	287	YieldSurface *YS_ ,
		288	PotentialSurface *PS_ ,
		289	EPState *EPS_,
		290	EvolutionLaw_S *ELS1_,
		291	EvolutionLaw_T *ELT1_ )
		292	:NDMaterial(tag, ND_TAG_Template3Dep)
		293	{
		294
		295	if ( YS_ )
		296	YS = YS_->newObj();
		297	else
		298	{
		299	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		300	exit(-1);
		301	}
		302
		303	if ( PS_ )
		304	PS = PS_->newObj();
		305	else
		306	{
		307	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		308	exit(-1);
		309	}
		310
		311	if ( EPS_ )
		312	EPS = EPS_->newObj();
		313	else
		314	{
		315	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		316	exit(-1);
		317	}
		318
		319	// Evolution laws
		320	if ( ELS1_ )
		321	ELS1 = ELS1_->newObj();
		322	else
		323	ELS1 = 0;
		324	ELS2 = 0;
		325	ELS3 = 0;
		326	ELS4 = 0;
		327
		328	if ( ELT1_ )
		329	ELT1 = ELT1_->newObj();
		330	else
		331	ELT1 = 0;
		332	ELT2 = 0;
		333	ELT3 = 0;
		334	ELT4 = 0;
		335	}
		336
		337	//================================================================================
		338	// Constructor 4
		339	// Two scalar evolution laws and one tensorial evolution law are provided!
		340	//================================================================================
		341
13	jeremic	342	Template3Dep::Template3Dep( int tag ,
12	jeremic	343	YieldSurface *YS_ ,
		344	PotentialSurface *PS_ ,
		345	EPState *EPS_,
		346	EvolutionLaw_S *ELS1_,
		347	EvolutionLaw_S *ELS2_,
		348	EvolutionLaw_T *ELT1_ )
		349	:NDMaterial(tag, ND_TAG_Template3Dep)
		350	{
		351
		352	if ( YS_ )
		353	YS = YS_->newObj();
		354	else
		355	{
		356	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		357	exit(-1);
		358	}
		359
		360	if ( PS_ )
		361	PS = PS_->newObj();
		362	else
		363	{
		364	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		365	exit(-1);
		366	}
		367
		368	if ( EPS_ )
		369	EPS = EPS_->newObj();
		370	else
		371	{
		372	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		373	exit(-1);
		374	}
		375
		376	if ( ELS1_ )
		377	ELS1 = ELS1_->newObj();
		378	else
		379	ELS1 = 0;
		380
		381	if ( ELS2_ )
		382	ELS2 = ELS2_->newObj();
		383	else
		384	ELS2 = 0;
		385
		386	ELS3 = 0;
		387	ELS4 = 0;
		388
		389	if ( ELT1_ )
		390	ELT1 = ELT1_->newObj();
		391	else
		392	ELT1 = 0;
13	jeremic	393
12	jeremic	394	ELT2 = 0;
		395	ELT3 = 0;
		396	ELT4 = 0;
		397	}
		398
		399	//================================================================================
		400	// Constructor 5
		401	// Two scalar evolution laws and two tensorial evolution law are provided!
		402	//================================================================================
		403
13	jeremic	404	Template3Dep::Template3Dep( int tag ,
12	jeremic	405	YieldSurface *YS_ ,
		406	PotentialSurface *PS_ ,
		407	EPState *EPS_,
		408	EvolutionLaw_S *ELS1_,
		409	EvolutionLaw_S *ELS2_,
		410	EvolutionLaw_T *ELT1_,
		411	EvolutionLaw_T *ELT2_)
		412	:NDMaterial(tag, ND_TAG_Template3Dep)
		413	{
		414
		415	if ( YS_ )
		416	YS = YS_->newObj();
		417	else
		418	{
		419	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		420	exit(-1);
		421	}
		422
		423	if ( PS_ )
		424	PS = PS_->newObj();
		425	else
		426	{
		427	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		428	exit(-1);
		429	}
		430
		431	if ( EPS_ )
		432	EPS = EPS_->newObj();
		433	else
		434	{
		435	g3ErrorHandler->fatal("Template3Dep:: Template3Dep failed to construct the template3Dep");
		436	exit(-1);
		437	}
		438
		439	if ( ELS1_ )
		440	ELS1 = ELS1_->newObj();
		441	else
		442	ELS1 = 0;
		443
		444	if ( ELS2_ )
		445	ELS2 = ELS2_->newObj();
		446	else
		447	ELS2 = 0;
		448	ELS3 = 0;
		449	ELS4 = 0;
		450
		451	if ( ELT1_ )
		452	ELT1 = ELT1_->newObj();
		453	else
		454	ELT1 = 0;
		455
		456	if ( ELT2_ )
		457	ELT2 = ELT2_->newObj();
		458	else
		459	ELT2 = 0;
		460	ELT3 = 0;
		461	ELT4 = 0;
		462	}
		463
		464	//================================================================================
		465	// Constructor 6
		466	// NO parameter is provided
		467	//================================================================================
		468
		469	Template3Dep::Template3Dep()
		470	:NDMaterial(0, ND_TAG_Template3Dep),ELS1(0), ELS2(0),ELS3(0), ELS4(0),
		471	ELT1(0), ELT2(0), ELT3(0), ELT4(0)
		472	{
		473	YS = new DPYieldSurface();
		474	PS = new DPPotentialSurface();
		475	EPS = new EPState();
		476	}
		477
		478
		479	//================================================================================
		480	// Destructor
		481	//================================================================================
		482
		483	Template3Dep::~Template3Dep()
365	jeremic	484	{
		485
130	jeremic	486	if (YS)
		487	delete YS;
12	jeremic	488
130	jeremic	489	if (PS)
		490	delete PS;
		491
		492	if (EPS)
		493	delete EPS;
		494
		495	if (ELS1)
		496	delete ELS1;
		497
		498	if (ELS2)
		499	delete ELS2;
		500
		501	if (ELS3)
		502	delete ELS3;
		503
		504	if (ELS4)
		505	delete ELS4;
		506
		507	if (ELT1)
		508	delete ELT1;
		509
		510	if (ELT2)
		511	delete ELT2;
		512
		513	if (ELT3)
		514	delete ELT3;
		515
		516	if (ELT4)
		517	delete ELT4;
365	jeremic	518
130	jeremic	519
12	jeremic	520	}
		521
		522	////================================================================================
		523	////copy constructor
		524	////================================================================================
		525	//Template3Dep::Template3Dep(const Template3Dep & rval) {
		526	//
		527	// YS = rval.YS->newObj();
		528	// PS = rval.PS->newObj();
		529	// EPS = rval.EPS->newObj();
		530	//
		531	// // Scalar Evolution Laws
		532	// if ( rval.getELS1() )
		533	// ELS1 = rval.getELS1()->newObj();
		534	// else
		535	// ELS1 = 0;
		536	//
		537	// if ( !rval.getELS2() )
		538	// ELS2 = 0;
		539	// else
		540	// ELS2 = rval.getELS2()->newObj();
		541	//
		542	// if ( !rval.getELS3() )
		543	// ELS3 = 0;
		544	// else
		545	// ELS3 = rval.getELS3()->newObj();
		546	//
		547	// if ( !rval.getELS4() )
		548	// ELS4 = 0;
		549	// else
		550	// ELS4 = rval.getELS4()->newObj();
		551	//
		552	// // Tensorial Evolution Laws
		553	// if ( rval.getELT1() )
		554	// ELT1 = rval.getELT1()->newObj();
		555	// else
		556	// ELT1 = 0;
		557	//
		558	// if ( !rval.getELT2() )
		559	// ELT2 = 0;
		560	// else
		561	// ELT2 = rval.getELT2()->newObj();
		562	//
		563	// if ( !rval.getELT3() )
		564	// ELT3 = 0;
		565	// else
		566	// ELT3 = rval.getELT3()->newObj();
		567	//
		568	// if ( !rval.getELT4() )
		569	// ELT4 = 0;
		570	// else
		571	// ELT4 = rval.getELT4()->newObj();
		572	//
		573	//}
		574	//
		575
		576
		577
		578	//================================================================================
		579	// Routine used to generate elastic compliance tensor D for this material point
		580	//================================================================================
		581	tensor Template3Dep::ElasticComplianceTensor(void) const
		582	{
		583	tensor ret( 4, def_dim_4, 0.0 );
		584
130	jeremic	585	double Ey = this->EPS->getEo();
		586	//cerr << " Eo= " << Ey;
12	jeremic	587	if (Ey == 0) {
		588	cout << endln << "Ey = 0! Can't give you D!!" << endln;
		589	exit(1);
		590	}
		591	double nuP =this->EPS->getnu();
		592
130	jeremic	593	//Update E
		594	stresstensor stc = (this->EPS)->getStress();
		595	double p = stc.p_hydrostatic();
		596	//cerr << " p = " << p;
		597
		598	double po = 100.0; //kPa
458	jeremic	599	if (p <= 0.001)
		600	p = 0.001;
365	jeremic	601	Ey = Ey * pow(p/po, 0.5); //0.5
130	jeremic	602	//cerr << " Ec = " << Ey << endln;
		603
12	jeremic	604	// Kronecker delta tensor
		605	tensor I2("I", 2, def_dim_2);
		606
		607	tensor I_ijkl = I2("ij")*I2("kl");
		608	//I_ijkl.null_indices();
		609	tensor I_ikjl = I_ijkl.transpose0110();
		610	tensor I_iljk = I_ijkl.transpose0111();
		611	tensor I4s = (I_ikjl+I_iljk)*0.5;
		612
		613	// Building compliance tensor
		614	ret = (I_ijkl * (-nuP/Ey)) + (I4s * ((1.0+nuP)/Ey));
		615
		616	return ret;
		617
		618	}
		619
		620	//================================================================================
		621	// Routine used to generate elastic stiffness tensor D for this material point
		622	//================================================================================
		623	tensor Template3Dep::ElasticStiffnessTensor(void) const
		624	{
		625	tensor ret( 4, def_dim_4, 0.0 );
130	jeremic	626
		627	double Ey = this->EPS->getEo();
		628	double nu =this->EPS->getnu();
		629
		630	//Update E
		631	stresstensor stc = (this->EPS)->getStress();
		632	double p = stc.p_hydrostatic();
		633	//cerr << " p = " << p;
12	jeremic	634
130	jeremic	635	double po = 100.0; //kPa
458	jeremic	636	if (p <= 0.001)
		637	p = 0.001;
365	jeremic	638	double E = Ey * pow(p/po, 0.5);//0.5
130	jeremic	639	//cerr << " Eo = " << Ey ;
		640	//cerr << " Ec = " << E << endln;
		641
13	jeremic	642
12	jeremic	643	// Kronecker delta tensor
		644	tensor I2("I", 2, def_dim_2);
		645
		646	tensor I_ijkl = I2("ij")*I2("kl");
		647
		648
		649	//I_ijkl.null_indices();
		650	tensor I_ikjl = I_ijkl.transpose0110();
		651	tensor I_iljk = I_ijkl.transpose0111();
		652	tensor I4s = (I_ikjl+I_iljk)*0.5;
		653
		654	//double x = I4s.trace();
		655	//cout << "xxxxxx " << x << endln;
		656
		657	//I4s.null_indices();
		658
		659	// Building elasticity tensor
		660	ret = I_ijkl( Enu / ( (1.0+nu)(1.0 - 2.0nu) ) ) + I4s*( E / (1.0 + nu) );
		661	//ret.null_indices();
		662
		663	return ret;
		664
		665
		666	}
		667
		668	//================================================================================
		669	int Template3Dep::setTrialStrain(const Vector &v)
		670	{
		671	// Not yet implemented
		672	return 0;
		673	}
		674
		675	//================================================================================
		676	int Template3Dep::setTrialStrain(const Vector &v, const Vector &r)
		677	{
		678	// Not yet implemented
		679	return 0;
		680	}
		681
		682	//================================================================================
		683	int Template3Dep::setTrialStrainIncr(const Vector &v)
		684	{
		685	// Not yet implemented
		686	return 0;
		687	}
		688
		689	//================================================================================
		690	int Template3Dep::setTrialStrainIncr(const Vector &v, const Vector &r)
		691	{
13	jeremic	692	//================================================================================
12	jeremic	693	// Not yet implemented
		694	return 0;
		695	}
		696
		697	//================================================================================
		698	const Matrix& Template3Dep::getTangent(void)
		699	{
		700	// Not yet implemented
		701	Matrix *M = new Matrix();
		702	return *M;
		703	}
		704
		705	//================================================================================
		706	const Vector& Template3Dep::getStress(void)
		707	{
		708	// Not yet implemented
		709	Vector *V = new Vector();
		710	return *V;
		711	}
		712
		713	//================================================================================
		714	const Vector& Template3Dep::getStrain(void)
		715	{
		716	// Not yet implemented
		717	Vector *V = new Vector();
		718	return *V;
		719	}
		720
		721	//================================================================================
130	jeremic	722	// what is the trial strain? Initial strain?
12	jeremic	723	int Template3Dep::setTrialStrain(const Tensor &v)
		724	{
		725	EPS->setStrain(v);
		726	return 0;
		727	}
		728
		729
		730	//================================================================================
		731	int Template3Dep::setTrialStrain(const Tensor &v, const Tensor &r)
		732	{
		733	EPS->setStrain(v);
		734	return 0;
		735	}
		736
		737	//================================================================================
		738
		739	int Template3Dep::setTrialStrainIncr(const Tensor &v)
		740	{
130	jeremic	741	//Need refining
		742
13	jeremic	743	EPState tmp_EPS = BackwardEulerEPState(v);
		744	if ( tmp_EPS.getConverged() ) {
		745	//setTrialEPS( tmp_EPS );
		746	setEPS( tmp_EPS );
130	jeremic	747	return 0;
13	jeremic	748	}
		749
130	jeremic	750	//int number_of_subincrements = 5;
		751	tmp_EPS = BESubIncrementation(v, NUM_OF_SUB_INCR);
13	jeremic	752	if ( tmp_EPS.getConverged() ) {
		753	//setTrialEPS( tmp_EPS );
		754	setEPS( tmp_EPS );
130	jeremic	755	return 0;
13	jeremic	756	}
		757	else {
130	jeremic	758	setEPS( tmp_EPS );
		759	return 1;
13	jeremic	760	}
130	jeremic	761
		762	//// for testing MD model only for no BE
		763	//EPState tmp_EPS = FESubIncrementation(v, NUM_OF_SUB_INCR);
		764	//setEPS( tmp_EPS );
		765	return 1;
12	jeremic	766	}
		767
		768	//================================================================================
		769	int Template3Dep::setTrialStrainIncr(const Tensor &v, const Tensor &r)
		770	{
		771	EPS->setStrain(v + EPS->getStrain() );
		772	return 0;
		773	}
		774
		775	//================================================================================
		776	const tensor& Template3Dep::getTangentTensor(void)
		777	{
13	jeremic	778	//Tensor Eep = EPS->getEep();
		779	return EPS->Eep;
12	jeremic	780	}
		781
		782	//================================================================================
146	jeremic	783	const stresstensor Template3Dep::getStressTensor(void)
12	jeremic	784	{
130	jeremic	785	//cout << *EPS;
		786	//stresstensor tmp;
		787	//tmp = EPS->getStress();
		788	//cout << "EPS->getStress() " << EPS->getStress() << endln;
		789
		790	//Something funny!!! happened here when returning EPS->getStress()
		791	// This function will return wrong Stress.
365	jeremic	792	stresstensor ret = EPS->getStress();
146	jeremic	793	return ret;
12	jeremic	794	}
		795
365	jeremic	796
12	jeremic	797	//================================================================================
365	jeremic	798	const straintensor Template3Dep::getStrainTensor(void)
12	jeremic	799	{
		800	return EPS->getStrain();
		801	}
		802
		803	//================================================================================
458	jeremic	804	const straintensor Template3Dep::getPlasticStrainTensor(void)
		805	{
		806	return EPS->getPlasticStrain();
		807	}
		808
		809
		810	//================================================================================
12	jeremic	811	int Template3Dep::commitState(void)
		812	{
458	jeremic	813	int err;
		814	err = getEPS()->commitState();
		815	return err;
12	jeremic	816	}
		817
		818	//================================================================================
		819	int Template3Dep::revertToLastCommit(void)
		820	{
130	jeremic	821	int err;
		822	err = EPS->revertToLastCommit();
		823	return err;
12	jeremic	824	}
		825	//================================================================================
		826	int Template3Dep::revertToStart(void)
		827	{
130	jeremic	828	int err;
		829	err = EPS->revertToStart();
		830	return err;
12	jeremic	831	}
		832
13	jeremic	833	////================================================================================
		834	//// Get one copy of the NDMaterial model
		835	////NDMaterial *Template3Dep::getCopy(void)
		836	//{
		837	// Template3Dep *t3d;
		838	// return t3d;
		839	//}
		840
12	jeremic	841	//================================================================================
13	jeremic	842	// Get one copy of the NDMaterial model
		843	NDMaterial * Template3Dep::getCopy(void)
12	jeremic	844	{
13	jeremic	845	NDMaterial * tmp =
12	jeremic	846	new Template3Dep( this->getTag() ,
		847	this->getYS() ,
		848	this->getPS() ,
		849	this->getEPS() ,
		850	this->getELS1() ,
		851	this->getELS2() ,
		852	this->getELS3() ,
		853	this->getELS4() ,
		854	this->getELT1() ,
		855	this->getELT2() ,
		856	this->getELT3() ,
		857	this->getELT4() );
		858
		859	return tmp;
		860
		861	}
		862
		863
		864	//================================================================================
13	jeremic	865	NDMaterial * Template3Dep::getCopy(const char *code)
12	jeremic	866	{
		867	if (strcmp(code,"Template3Dep") == 0)
		868	{
133	fmk	869	Template3Dep * tmp =
12	jeremic	870	new Template3Dep( this->getTag() ,
		871	this->getYS() ,
		872	this->getPS() ,
		873	this->getEPS() ,
		874	this->getELS1() ,
		875	this->getELS2() ,
		876	this->getELS3() ,
		877	this->getELS4() ,
		878	this->getELT1() ,
		879	this->getELT2() ,
		880	this->getELT3() ,
		881	this->getELT4() );
		882
13	jeremic	883	tmp->getEPS()->setInit();
12	jeremic	884	return tmp;
		885	}
		886	else
		887	{
		888	g3ErrorHandler->fatal("Template3Dep::getCopy failed to get model %s", code);
		889	return 0;
		890	}
		891
		892	}
		893
		894	//================================================================================
		895	const char *Template3Dep::getType(void) const
		896	{
		897	return "Template3Dep";
		898	}
		899
		900	//================================================================================
130	jeremic	901	//??What is the Order????????? might be the
12	jeremic	902
		903	int Template3Dep::getOrder(void) const
		904	{
130	jeremic	905	return 6;
12	jeremic	906	}
		907
		908	//================================================================================
		909	int Template3Dep::sendSelf(int commitTag, Channel &theChannel)
		910	{
		911	// Not yet implemented
		912	return 0;
		913	}
		914
		915	//================================================================================
		916	int Template3Dep::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
		917	{
		918	// Not yet implemented
		919	return 0;
		920	}
		921
		922	//================================================================================
		923	void
		924	Template3Dep::Print(ostream &s, int flag)
		925	{
		926	s << (*this);
		927	}
		928
		929
		930
		931	//private utilities
		932
		933	//================================================================================
		934	// Set new EPState
		935	//================================================================================
130	jeremic	936
		937	void Template3Dep::setEPS( EPState & rval)
		938	{
		939	//EPState eps = rval; older buggy one
		940	//EPS = rval.newObj();
		941	/*
		942	//EPS->setEo(rval.getEo());
		943	EPS->setE(rval.getE());
		944	//EPS->setnu(getnu());
		945	//EPS->setrho(getrho());
		946	EPS->setStress(rval.getStress());
		947	EPS->setStrain(rval.getStrain());
		948	EPS->setElasticStrain(rval.getElasticStrain());
		949	EPS->setPlasticStrain(rval.getPlasticStrain());
		950	EPS->setdElasticStrain(rval.getdElasticStrain());
		951	EPS->setdPlasticStrain(rval.getdPlasticStrain());
		952	//EPS->setNScalarVar(rval.getNScalarVar());
		953	for (int i = 0; i <rval.getNScalarVar(); i++)
		954	EPS->setScalarVar(i, rval.getScalarVar(i));
		955
		956
		957	//EPS->setNTensorVar(rval.getNTensorVar());
		958	EPS->setTensorVar(rval.getTensorVar());
		959	EPS->setEep(rval.getEep());
		960	EPS->Stress_commit=(rval.getStress_commit());
		961	EPS->Strain_commit=(rval.getStrain_commit());
		962
		963	for (int i = 0; i <rval.getNScalarVar(); i++)
		964	EPS->ScalarVar_commit[i] = rval.getScalarVar_commit()[i];
		965
		966	for (int i = 0; i <rval.getNTensorVar(); i++)
		967	EPS->TensorVar_commit[i] = (rval.getTensorVar_commit()[i]);
		968
		969	EPS->Eep_commit = (rval.getEep_commit());
		970	//EPS->setStress_init(getStress_init());
		971	//EPS->setStrain_init(getStrain_init());
		972	//EPS->setScalarVar_init(getScalarVar_init());
		973	//EPS->setTensorVar_init(getTensorVar_init());
		974	//EPS->setEep_init(getEep_init());
		975	EPS->setConverged(rval.getConverged());
		976
		977	*/
		978
		979	EPS->setEo(rval.getEo());
		980	EPS->setE(rval.getE());
		981	EPS->setStress(rval.getStress());
		982	EPS->setStrain(rval.getStrain());
		983	EPS->setElasticStrain(rval.getElasticStrain());
		984	EPS->setPlasticStrain(rval.getPlasticStrain());
		985	EPS->setdElasticStrain(rval.getdElasticStrain());
		986	EPS->setdPlasticStrain(rval.getdPlasticStrain());
		987	EPS->setNScalarVar( rval.getNScalarVar() );
		988
458	jeremic	989	for (int i = 0; i <rval.getNScalarVar(); i++)
130	jeremic	990	EPS->setScalarVar(i+1, rval.getScalarVar(i+1));
		991
		992
		993	EPS->setNTensorVar( rval.getNTensorVar() );
		994	EPS->setTensorVar(rval.getTensorVar());
		995	EPS->setEep(rval.getEep());
		996	EPS->setStress_commit(rval.getStress_commit());
		997	EPS->setStrain_commit(rval.getStrain_commit());
		998
458	jeremic	999	for (int i = 0; i <rval.getNScalarVar(); i++)
130	jeremic	1000	EPS->setScalarVar_commit(i+1, rval.getScalarVar_commit(i+1));
		1001
458	jeremic	1002	for (int i = 0; i <rval.getNTensorVar(); i++)
130	jeremic	1003	EPS->setTensorVar_commit(i+1, rval.getTensorVar_commit(i+1));
		1004
		1005	EPS->Eep_commit = (rval.getEep_commit());
		1006	EPS->Stress_init = rval.getStress_init();
		1007	EPS->Strain_init = rval.getStrain_init();
		1008
458	jeremic	1009	for (int i = 0; i <rval.getNScalarVar(); i++)
130	jeremic	1010	EPS->setScalarVar_init(i+1, rval.getScalarVar_init(i+1));
		1011
458	jeremic	1012	for (int i = 0; i <rval.getNTensorVar(); i++)
130	jeremic	1013	EPS->setTensorVar_init(i+1, rval.getTensorVar_init(i+1));
		1014
		1015	EPS->Eep_init = rval.getEep_init();
		1016	EPS->setConverged(rval.getConverged());
		1017
		1018
12	jeremic	1019	}
		1020
13	jeremic	1021
130	jeremic	1022
12	jeremic	1023	//================================================================================
		1024	// Get the Yield Surface
		1025	//================================================================================
		1026	YieldSurface * Template3Dep::getYS() const
		1027	{
		1028	return YS;
		1029	}
		1030
		1031
		1032	//================================================================================
		1033	// Get the Potential Surface
		1034	//================================================================================
		1035	PotentialSurface * Template3Dep::getPS() const
		1036	{
		1037	return PS;
		1038	}
		1039
		1040	//================================================================================
		1041	// Get the EPState
		1042	//================================================================================
		1043	EPState * Template3Dep::getEPS() const
		1044	{
		1045	return EPS;
		1046	}
		1047
13	jeremic	1048
12	jeremic	1049	//================================================================================
		1050	// Get the 1st Salar evolution law
		1051	//================================================================================
		1052
		1053	EvolutionLaw_S * Template3Dep::getELS1() const
		1054	{
		1055	return ELS1;
		1056	}
		1057
		1058	//================================================================================
		1059	// Get the 2ndst Salar evolution law
		1060	//================================================================================
		1061	EvolutionLaw_S * Template3Dep::getELS2() const
		1062	{
		1063	return ELS2;
		1064	}
		1065
		1066	//================================================================================
		1067	// Get the 2ndst Salar evolution law
		1068	//================================================================================
		1069	EvolutionLaw_S * Template3Dep::getELS3() const
		1070	{
		1071	return ELS3;
		1072	}
		1073	//================================================================================
		1074	// Get the 2ndst Salar evolution law
		1075	//================================================================================
		1076	EvolutionLaw_S * Template3Dep::getELS4() const
		1077	{
		1078	return ELS4;
		1079	}
		1080
		1081
		1082	//================================================================================
		1083	// Get the 1st tensorial evolution law
		1084	//================================================================================
		1085
		1086	EvolutionLaw_T * Template3Dep::getELT1() const
		1087	{
		1088	return ELT1;
		1089	}
		1090
		1091	//================================================================================
		1092	// Get the 2nd tensorial evolution law
		1093	//================================================================================
		1094
		1095	EvolutionLaw_T * Template3Dep::getELT2() const
		1096	{
		1097	return ELT2;
		1098	}
		1099	//================================================================================
		1100	// Get the 3rd tensorial evolution law
		1101	//================================================================================
		1102
		1103	EvolutionLaw_T * Template3Dep::getELT3() const
		1104	{
		1105	return ELT3;
		1106	}
		1107	//================================================================================
		1108	// Get the 4th tensorial evolution law
		1109	//================================================================================
		1110
		1111	EvolutionLaw_T * Template3Dep::getELT4() const
		1112	{
		1113	return ELT4;
		1114	}
		1115
		1116
		1117	//================================================================================
13	jeremic	1118	// Predictor EPState by Forward, Backward, MidPoint Methods...
		1119	//================================================================================
		1120
		1121	EPState Template3Dep::PredictorEPState(straintensor & strain_increment)
		1122	{
		1123
		1124	EPState Predictor = ForwardEulerEPState( strain_increment);
		1125	//EPState Predictor = SemiBackwardEulerEPState( strain_increment, material_point);
		1126	return Predictor;
		1127
		1128	}
		1129
		1130	//================================================================================
		1131	// New EPState using Forward Euler Algorithm
		1132	//================================================================================
		1133	EPState Template3Dep::ForwardEulerEPState( straintensor &strain_increment)
		1134	{
130	jeremic	1135	// Volumetric strain
		1136	double st_vol = strain_increment.p_hydrostatic();
		1137
13	jeremic	1138	//EPState forwardEPS( *(material_point->getEPS()) );
		1139	EPState forwardEPS( *(this->getEPS()) );
		1140	//cout <<"start eps: " << forwardEPS;
365	jeremic	1141	//cout << "\nForwardEulerEPState strain_increment " << strain_increment << endln;
13	jeremic	1142
		1143	// Building elasticity tensor
		1144	tensor E = ElasticStiffnessTensor();
130	jeremic	1145	//tensor Eep = ElasticStiffnessTensor();
		1146	tensor Eep = E;
13	jeremic	1147	tensor D = ElasticComplianceTensor();
		1148	E.null_indices();
		1149	D.null_indices();
		1150
		1151	//Checking E and D
		1152	//tensor ed = E("ijpq") * D("pqkl");
		1153	//double edd = ed.trace(); // = 3.
		1154
		1155	stresstensor stress_increment = E("ijpq") * strain_increment("pq");
		1156	stress_increment.null_indices();
		1157	//cout << " stress_increment: " << stress_increment << endln;
		1158
		1159	EPState startEPS( *(getEPS()) );
		1160	stresstensor start_stress = startEPS.getStress();
		1161	start_stress.null_indices();
		1162	//cout << "===== start_EPS =====: " << startEPS;
		1163
		1164	double f_start = 0.0;
		1165	double f_pred = 0.0;
		1166
		1167	EPState IntersectionEPS( startEPS );
		1168
		1169	EPState ElasticPredictorEPS( startEPS );
		1170	stresstensor elastic_predictor_stress = start_stress + stress_increment;
		1171	ElasticPredictorEPS.setStress( elastic_predictor_stress );
		1172	//cout << " Elastic_predictor_stress: " << elastic_predictor_stress << endln;
		1173
		1174	f_start = this->getYS()->f( &startEPS );
		1175	//::printf("\n############## f_start = %.10e ",f_start);
365	jeremic	1176	//cout << "\n####### f_start = " << f_start;
13	jeremic	1177
		1178	f_pred = this->getYS()->f( &ElasticPredictorEPS );
		1179	//::printf("############## f_pred = %.10e\n\n",f_pred);
365	jeremic	1180	//cout << " ####### f_pred = " << f_pred << "\n";
13	jeremic	1181
		1182	stresstensor intersection_stress = start_stress; // added 20april2000 for forward euler
		1183	stresstensor elpl_start_stress = start_stress;
		1184	stresstensor true_stress_increment = stress_increment;
		1185
130	jeremic	1186	if ( f_start <= 0 && f_pred <= 0 \|\| f_start > f_pred )
13	jeremic	1187	{
		1188	//Updating elastic strain increment
		1189	straintensor estrain = ElasticPredictorEPS.getElasticStrain();
		1190	straintensor tstrain = ElasticPredictorEPS.getStrain();
		1191	estrain = estrain + strain_increment;
		1192	tstrain = tstrain + strain_increment;
		1193	ElasticPredictorEPS.setElasticStrain( estrain );
		1194	ElasticPredictorEPS.setStrain( tstrain );
		1195	ElasticPredictorEPS.setdElasticStrain( strain_increment );
130	jeremic	1196
		1197	//Evolve parameters like void ratio (e) according to elastic strain and elastic stress--for MD model especially
13	jeremic	1198	//double Delta_lambda = 0.0;
		1199	//material_point.EL->UpdateVar( &ElasticPredictorEPS, 1);
130	jeremic	1200	// Update E_Young and e according to current stress state before evaluate ElasticStiffnessTensor
		1201	if ( getELT1() )
		1202	getELT1()->updateEeDm(&ElasticPredictorEPS, st_vol, 0.0);
13	jeremic	1203
		1204	//cout <<" strain_increment.Iinvariant1() " << strain_increment.Iinvariant1() << endln;
		1205
		1206	ElasticPredictorEPS.setEep(E);
		1207	return ElasticPredictorEPS;
		1208	}
		1209
		1210	if ( f_start <= 0 && f_pred > 0 )
		1211	{
		1212	intersection_stress =
		1213	yield_surface_cross( start_stress, elastic_predictor_stress);
365	jeremic	1214	//cout << " start_stress: " << start_stress << endln;
13	jeremic	1215	//cout << " Intersection_stress: " << intersection_stress << endln;
		1216
		1217	IntersectionEPS.setStress( intersection_stress );
		1218	//intersection_stress.reportshort("Intersection stress \n");
		1219
		1220	elpl_start_stress = intersection_stress;
		1221	//elpl_start_stress.reportshortpqtheta("elpl start stress AFTER \n");
		1222
		1223	true_stress_increment = elastic_predictor_stress - elpl_start_stress;
		1224	//true_stress_increment.null_indices();
130	jeremic	1225
		1226	if ( getELT1() )
		1227	getELT1()->updateEeDm(&ElasticPredictorEPS, st_vol, 0.0);
13	jeremic	1228
		1229	}
		1230
		1231
		1232	forwardEPS.setStress( elpl_start_stress );
		1233
		1234	//cout <<"elpl start eps: " << forwardEPS;
130	jeremic	1235	//double f_cross = this->getYS()->f( &forwardEPS );
13	jeremic	1236	//cout << " ####### f_cross = " << f_cross << "\n";
		1237
		1238	//set the initial value of D once the current stress hits the y.s. for Manzari-Dafalias Model
		1239	//if ( f_start <= 0 && f_pred > 0 )
		1240	// material_point.EL->setInitD(&forwardEPS);
		1241	//cout << " inside ConstitutiveDriver after setInitD " << forwardEPS;
		1242
		1243
		1244	// pulling out some tensor and double definitions
		1245	//tensor dFods( 2, def_dim_2, 0.0);
		1246	//tensor dQods( 2, def_dim_2, 0.0);
		1247	stresstensor dFods;
		1248	stresstensor dQods;
		1249	// stresstensor s; // deviator
		1250	tensor H( 2, def_dim_2, 0.0);
		1251	tensor temp1( 2, def_dim_2, 0.0);
		1252	tensor temp2( 2, def_dim_2, 0.0);
		1253	double lower = 0.0;
		1254	tensor temp3( 2, def_dim_2, 0.0);
		1255
		1256	double Delta_lambda = 0.0;
		1257	double h_s[4] = {0.0, 0.0, 0.0, 0.0};
		1258	double xi_s[4] = {0.0, 0.0, 0.0, 0.0};
		1259	stresstensor h_t[4];
		1260	stresstensor xi_t[4];
		1261	double hardMod_ = 0.0;
		1262
		1263	//double Dq_ast = 0.0;
		1264	//double q_ast_entry = 0.0;
		1265	//double q_ast = 0.0;
		1266
		1267	stresstensor plastic_stress;
		1268	straintensor plastic_strain;
		1269	stresstensor elastic_plastic_stress;
		1270	// ::printf("\n��...... felpred = %lf\n",felpred);
		1271
		1272	if ( f_pred >= 0 ) {
		1273
		1274
		1275	dFods = getYS()->dFods( &forwardEPS );
		1276	dQods = getPS()->dQods( &forwardEPS );
		1277
		1278	//cout << "dF/ds" << dFods << endln;
		1279	//cout << "dQ/ds" << dQods << endln;
		1280
		1281	// Tensor H_kl ( eq. 5.209 ) W.F. Chen
		1282	H = E("ijkl")dQods("kl"); //E_ijkl R_kl
		1283	H.null_indices();
		1284	temp1 = dFods("ij") * E("ijkl"); // L_ij * E_ijkl
		1285	temp1.null_indices();
		1286	temp2 = temp1("ij")dQods("ij"); // L_ij E_ijkl * R_kl
		1287	temp2.null_indices();
		1288	lower = temp2.trace();
		1289
		1290	// Evaluating the hardening modulus: sum of (df/dq) qbar
		1291
		1292	hardMod_ = 0.0;
		1293	//Of 1st scalar internal vars
		1294	if ( getELS1() ) {
		1295	h_s[0] = getELS1()->h_s(&forwardEPS, getPS());
		1296	xi_s[0] = getYS()->xi_s1( &forwardEPS );
		1297	hardMod_ = hardMod_ + h_s[0] * xi_s[0];
		1298	}
		1299
		1300	//Of 2nd scalar internal vars
		1301	if ( getELS2() ) {
		1302	h_s[1] = getELS2()->h_s( &forwardEPS, getPS());
		1303	xi_s[1] = getYS()->xi_s2( &forwardEPS );
		1304	hardMod_ = hardMod_ + h_s[1] * xi_s[1];
		1305	}
		1306
		1307	//Of 3rd scalar internal vars
		1308	if ( getELS3() ) {
		1309	h_s[2] = getELS3()->h_s( &forwardEPS, getPS());
		1310	xi_s[2] = getYS()->xi_s3( &forwardEPS );
		1311	hardMod_ = hardMod_ + h_s[2] * xi_s[2];
		1312	}
		1313
		1314	//Of 4th scalar internal vars
		1315	if ( getELS4() ) {
		1316	h_s[3] = getELS4()->h_s( &forwardEPS, getPS());
		1317	xi_s[3] = getYS()->xi_s4( &forwardEPS );
		1318	hardMod_ = hardMod_ + h_s[3] * xi_s[3];
		1319	}
		1320
		1321	//Of tensorial internal var
		1322	// 1st tensorial var
		1323	if ( getELT1() ) {
		1324	h_t[0] = getELT1()->h_t(&forwardEPS, getPS());
		1325	xi_t[0] = getYS()->xi_t1( &forwardEPS );
		1326	tensor hm = (h_t[0])("ij") * (xi_t[0])("ij");
		1327	hardMod_ = hardMod_ + hm.trace();
		1328	}
		1329
		1330	// 2nd tensorial var
		1331	if ( getELT2() ) {
		1332	h_t[1] = getELT2()->h_t( &forwardEPS, getPS());
		1333	xi_t[1] = getYS()->xi_t2( &forwardEPS );
		1334	tensor hm = (h_t[1])("ij") * (xi_t[1])("ij");
		1335	hardMod_ = hardMod_ + hm.trace();
		1336	}
		1337
		1338	// 3rd tensorial var
		1339	if ( getELT3() ) {
		1340	h_t[2] = getELT3()->h_t( &forwardEPS, getPS());
		1341	xi_t[2] = getYS()->xi_t3( &forwardEPS );
		1342	tensor hm = (h_t[2])("ij") * (xi_t[2])("ij");
		1343	hardMod_ = hardMod_ + hm.trace();
		1344	}
		1345
		1346	// 4th tensorial var
		1347	if ( getELT4() ) {
		1348	h_t[3] = getELT4()->h_t(&forwardEPS, getPS());
		1349	xi_t[3] = getYS()->xi_t4( &forwardEPS );
		1350	tensor hm = (h_t[3])("ij") * (xi_t[3])("ij");
		1351	hardMod_ = hardMod_ + hm.trace();
		1352	}
		1353
		1354	// Subtract accumulated hardMod_ from lower
		1355	lower = lower - hardMod_;
		1356
		1357	//Calculating Kp according to Kp = - (df/dq) qbar
		1358	//double Kp = material_point.EL->getKp(&forwardEPS, norm_dQods);
		1359	//Kp = 0.0;
		1360	//cout << endln << ">>>>>>>>> Lower = " << lower << endln;
		1361	//lower = lower + Kp;
		1362	//cout << endln << ">>>>>>>>> Kp = " << Kp << endln;
		1363
		1364	//cout << " stress_increment "<< stress_increment << endln;
		1365	//cout << " true_stress_increment "<< true_stress_increment << endln;
		1366
		1367	temp3 = dFods("ij") * true_stress_increment("ij"); // L_ij * E_ijkl * d e_kl (true ep strain increment)
		1368	temp3.null_indices();
		1369	//cout << " temp3.trace() " << temp3.trace() << endln;
		1370	Delta_lambda = (temp3.trace())/lower;
		1371	//cout << "Delta_lambda " << Delta_lambda << endln;
		1372	//if (Delta_lambda<0.0) Delta_lambda=0.0;
		1373
		1374	plastic_stress = H("kl") * Delta_lambda;
		1375	plastic_strain = dQods("kl") * Delta_lambda; // plastic strain increment
		1376	plastic_stress.null_indices();
		1377	plastic_strain.null_indices();
		1378	//cout << "plastic_stress = " << plastic_stress << endln;
		1379	//cout << "plastic_strain = " << plastic_strain << endln;
		1380	//cout << "plastic_strain I1= " << plastic_strain.Iinvariant1() << endln;
		1381	//cout << "plastic_strain vol " << Delta_lambda * ( forwardEPS.getScalarVar( 2 ) )<< endln ;
		1382	//cout << " q=" << Delta_lambda * dQods.q_deviatoric()<< endln;
		1383	//plastic_stress.reportshort("plastic stress (with delta_lambda)\n");
		1384
		1385	elastic_plastic_stress = elastic_predictor_stress - plastic_stress;
		1386	//elastic_plastic_stress.reportshortpqtheta("FE elastic plastic stress \n");
		1387
		1388	//calculating elatic strain increment
		1389	//stresstensor dstress_el = elastic_plastic_stress - start_stress;
		1390	//straintensor elastic_strain = D("ijpq") * dstress_el("pq");
		1391	straintensor elastic_strain = strain_increment - plastic_strain; // elastic strain increment
		1392	//cout << "elastic_strain I1=" << elastic_strain.Iinvariant1() << endln;
		1393	//cout << "elastic_strain " << elastic_strain << endln;
		1394	//cout << "strain increment I1=" << strain_increment.Iinvariant1() << endln;
		1395	//cout << "strain increment " << strain_increment << endln;
		1396
		1397	straintensor estrain = forwardEPS.getElasticStrain(); //get old elastic strain
		1398	straintensor pstrain = forwardEPS.getPlasticStrain(); //get old plastic strain
		1399
		1400	straintensor tstrain = forwardEPS.getStrain(); //get old total strain
		1401	pstrain = pstrain + plastic_strain;
		1402	estrain = estrain + elastic_strain;
		1403	tstrain = tstrain + elastic_strain + plastic_strain;
		1404
		1405	//Setting de_p, de_e, total plastic, elastic strain, and total strain
		1406	forwardEPS.setdPlasticStrain( plastic_strain );
		1407	forwardEPS.setdElasticStrain( elastic_strain );
		1408	forwardEPS.setPlasticStrain( pstrain );
		1409	forwardEPS.setElasticStrain( estrain );
		1410	forwardEPS.setStrain( tstrain );
		1411
		1412	//================================================================
		1413	//Generating Eep using dQods at the intersection point
		1414	dFods = getYS()->dFods( &IntersectionEPS );
		1415	dQods = getPS()->dQods( &IntersectionEPS );
		1416
		1417	tensor upperE1 = E("pqkl")*dQods("kl");
		1418	upperE1.null_indices();
		1419	tensor upperE2 = dFods("ij")*E("ijmn");
		1420	upperE2.null_indices();
		1421	tensor upperE = upperE1("pq") * upperE1("mn");
		1422	upperE.null_indices();
		1423
		1424	temp2 = upperE2("ij")dQods("ij"); // L_ij E_ijkl * R_kl
		1425	temp2.null_indices();
		1426	lower = temp2.trace();
		1427
		1428	// Evaluating the hardening modulus: sum of (df/dq) qbar
		1429
		1430	hardMod_ = 0.0;
		1431	//Of 1st scalar internal vars
		1432	if ( getELS1() ) {
		1433	h_s[0] = getELS1()->h_s( &IntersectionEPS, getPS());
		1434	xi_s[0] = getYS()->xi_s1( &IntersectionEPS );
		1435	hardMod_ = hardMod_ + h_s[0] * xi_s[0];
		1436	}
		1437
		1438	//Of 2nd scalar internal vars
		1439	if ( getELS2() ) {
		1440	h_s[1] = getELS2()->h_s( &IntersectionEPS, getPS());
		1441	xi_s[1] = getYS()->xi_s2( &IntersectionEPS );
		1442	hardMod_ = hardMod_ + h_s[1] * xi_s[1];
		1443	}
		1444
		1445	//Of 3rd scalar internal vars
		1446	if ( getELS3() ) {
		1447	h_s[2] = getELS3()->h_s( &IntersectionEPS, getPS());
		1448	xi_s[2] = getYS()->xi_s3( &IntersectionEPS );
		1449	hardMod_ = hardMod_ + h_s[2] * xi_s[2];
		1450	}
		1451
		1452	//Of 4th scalar internal vars
		1453	if ( getELS4() ) {
		1454	h_s[3] = getELS4()->h_s( &IntersectionEPS, getPS());
		1455	xi_s[3] = getYS()->xi_s4( &IntersectionEPS );
		1456	hardMod_ = hardMod_ + h_s[3] * xi_s[3];
		1457	}
		1458
		1459	//Of tensorial internal var
		1460	// 1st tensorial var
		1461	if ( getELT1() ) {
		1462	h_t[0] = getELT1()->h_t(&IntersectionEPS, getPS());
		1463	xi_t[0] = getYS()->xi_t1( &IntersectionEPS );
		1464	tensor hm = (h_t[0])("ij") * (xi_t[0])("ij");
		1465	hardMod_ = hardMod_ + hm.trace();
		1466	}
		1467
		1468	// 2nd tensorial var
		1469	if ( getELT2() ) {
		1470	h_t[1] = getELT2()->h_t( &IntersectionEPS, getPS());
		1471	xi_t[1] = getYS()->xi_t2( &IntersectionEPS );
		1472	tensor hm = (h_t[1])("ij") * (xi_t[1])("ij");
		1473	hardMod_ = hardMod_ + hm.trace();
		1474	}
		1475
		1476	// 3rd tensorial var
		1477	if ( getELT3() ) {
		1478	h_t[2] = getELT3()->h_t(&IntersectionEPS, getPS());
		1479	xi_t[2] = getYS()->xi_t3( &IntersectionEPS );
		1480	tensor hm = (h_t[2])("ij") * (xi_t[2])("ij");
		1481	hardMod_ = hardMod_ + hm.trace();
		1482	}
		1483
		1484	// 4th tensorial var
		1485	if ( getELT4() ) {
		1486	h_t[3] = getELT4()->h_t( &IntersectionEPS, getPS());
		1487	xi_t[3] = getYS()->xi_t4( &IntersectionEPS );
		1488	tensor hm = (h_t[3])("ij") * (xi_t[3])("ij");
		1489	hardMod_ = hardMod_ + hm.trace();
		1490	}
		1491
		1492	// Subtract accumulated hardMod_ from lower
		1493	lower = lower - hardMod_;
		1494
		1495	tensor Ep = upperE*(1./lower);
		1496
		1497	// elastoplastic constitutive tensor
		1498	Eep = Eep - Ep;
		1499
		1500	//cout <<" after calculation---Eep.rank()= " << Eep.rank() <<endln;
130	jeremic	1501	//Eep.printshort(" IN template ");
13	jeremic	1502
		1503	//--// before the surface is been updated !
		1504	//--// f_Final = Criterion.f(elastic_plastic_stress);
		1505	//--
		1506	//-- q_ast_entry = Criterion.kappa_get(elastic_plastic_stress);
		1507	//--
		1508	//--// h_ = h(elastic_plastic_stress);
		1509	//-- Dq_ast = Delta_lambda * h_ * just_this_PP;
		1510	//--// if (Dq_ast < 0.0 ) Dq_ast = 0.0;
		1511	//--// Dq_ast = fabs(Delta_lambda * h_ * just_this_PP); // because of softening response...
		1512	//--//::printf(" h_=%.6e q_ast_entry=%.6e Dq_ast=%.6e Delta_lambda=%.6e\n", h_, q_ast_entry, Dq_ast, Delta_lambda);
		1513	//--
		1514	//-- current_lambda_set(Delta_lambda);
		1515	//--
		1516	//-- q_ast = q_ast_entry + Dq_ast;
		1517	//-- Criterion.kappa_set( elastic_plastic_stress, q_ast);
		1518	//--//::fprintf(stdout," Criterion.kappa_get(elastic_plastic_stress)=%.8e\n",Criterion.kappa_get(elastic_plastic_stress));
		1519	//--//::fprintf(stderr," Criterion.kappa_get(elastic_plastic_stress)=%.8e\n",Criterion.kappa_get(elastic_plastic_stress));
		1520	//--
		1521	//--
		1522	//--//::fprintf(stdout," ######## predictor --> q_ast_entry=%.8e Dq_ast=%.8e q_ast=%.8e\n",q_ast_entry, Dq_ast, q_ast);
		1523	//--//::fprintf(stderr," ######## predictor --> q_ast_entry=%.8e Dq_ast=%.8e q_ast=%.8e\n",q_ast_entry, Dq_ast, q_ast);
		1524	//--
		1525	//--//::fprintf(stdout,"ForwardEulerStress IN Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress));
		1526	//--//::fprintf(stderr,"ForwardEulerStress IN Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress));
		1527	//--
		1528
		1529	//new EPState in terms of stress
		1530	forwardEPS.setStress(elastic_plastic_stress);
		1531	//cout <<"inside constitutive driver: forwardEPS "<< forwardEPS;
		1532
		1533	forwardEPS.setEep(Eep);
		1534	//forwardEPS.getEep().printshort(" after set");
		1535
		1536	//Before update all the internal vars
		1537	double f_forward = getYS()->f( &forwardEPS );
		1538	//::printf("\n************ Before f_forward = %.10e\n",f_forward);
		1539
		1540	//Evolve the surfaces and hardening vars
		1541	int NS = forwardEPS.getNScalarVar();
		1542	int NT = forwardEPS.getNTensorVar();
		1543
		1544	double dS= 0.0;
		1545	double S = 0.0;
		1546
458	jeremic	1547	for (int ii = 1; ii <= NS; ii++) {
13	jeremic	1548	dS = Delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var
		1549	S = forwardEPS.getScalarVar(ii); // Get the old value of the scalar internal var
		1550	forwardEPS.setScalarVar(ii, S + dS ); // Update internal scalar var
		1551	}
		1552
		1553	stresstensor dT;
		1554	stresstensor T;
		1555	stresstensor new_T;
		1556
458	jeremic	1557	for (int ii = 1; ii <= NT; ii++) {
157	jeremic	1558	dT = h_t[ii-1]*Delta_lambda ; // Increment to the tensor internal var
13	jeremic	1559	T = forwardEPS.getTensorVar(ii); // Get the old value of the tensor internal var
		1560	new_T = T + dT;
		1561	forwardEPS.setTensorVar(ii, new_T );
		1562	}
		1563
130	jeremic	1564	// Update E_Young and e according to current stress state before evaluate ElasticStiffnessTensor
		1565	int err = 0;
		1566	if ( getELT1() )
		1567	err = getELT1()->updateEeDm(&forwardEPS, st_vol, Delta_lambda);
		1568
13	jeremic	1569	//tensor tempx = plastic_strain("ij") * plastic_strain("ij");
		1570	//double tempxd = tempx.trace();
		1571	//double e_eq = pow( 2.0 * tempxd / 3.0, 0.5 );
		1572	////cout << "e_eq = " << e_eq << endln;
		1573	//
		1574	//double dalfa1 = e_eq * 10;
		1575	//double alfa1 = forwardEPS.getScalarVar(1);
		1576
		1577
		1578
		1579	//cout << "UpdateAllVars " << forwardEPS<< endln;
		1580
		1581	//After update all the internal vars
		1582	f_forward = getYS()->f( &forwardEPS );
		1583	//::printf("\n************ After f_forward = %.10e\n\n",f_forward);
		1584
		1585
		1586	}
		1587
		1588	//::fprintf(stderr,"ForwardEulerStress EXIT Criterion.kappa_get(start_stress)=%.8e\n",Criterion.kappa_get(start_stress));
		1589	return forwardEPS;
		1590	}
		1591
		1592	//================================================================================
		1593	// Starting EPState using Semi Backward Euler Starting Point
		1594	//================================================================================
		1595	EPState Template3Dep::SemiBackwardEulerEPState( const straintensor &strain_increment)
		1596	{
		1597	stresstensor start_stress;
		1598	EPState SemibackwardEPS( *(this->getEPS()) );
		1599	start_stress = SemibackwardEPS.getStress();
		1600
		1601	// building elasticity tensor
		1602	//tensor E = Criterion.StiffnessTensorE(Ey,nu);
		1603	tensor E = ElasticStiffnessTensor();
		1604	// building compliance tensor
		1605	// tensor D = Criterion.ComplianceTensorD(Ey,nu);
		1606
		1607	// pulling out some tensor and double definitions
		1608	tensor dFods(2, def_dim_2, 0.0);
		1609	tensor dQods(2, def_dim_2, 0.0);
		1610	tensor temp2(2, def_dim_2, 0.0);
		1611	double lower = 0.0;
		1612	double Delta_lambda = 0.0;
		1613
		1614	EPState E_Pred_EPS( *(this->getEPS()) );
		1615
		1616	straintensor strain_incr = strain_increment;
		1617	stresstensor stress_increment = E("ijkl")*strain_incr("kl");
		1618	stress_increment.null_indices();
		1619	// stress_increment.reportshort("stress Increment\n");
		1620
		1621
		1622	stresstensor plastic_stress;
		1623	stresstensor elastic_predictor_stress;
		1624	stresstensor elastic_plastic_stress;
		1625	//.. double Felplpredictor = 0.0;
		1626
		1627	double h_s[4] = {0.0, 0.0, 0.0, 0.0};
		1628	double xi_s[4] = {0.0, 0.0, 0.0, 0.0};
		1629	stresstensor h_t[4];
		1630	stresstensor xi_t[4];
		1631	double hardMod_ = 0.0;
		1632
		1633	double S = 0.0;
		1634	double dS = 0.0;
		1635	stresstensor T;
		1636	stresstensor dT;
		1637	//double Dq_ast = 0.0;
		1638	//double q_ast_entry = 0.0;
		1639	//double q_ast = 0.0;
		1640
		1641	elastic_predictor_stress = start_stress + stress_increment;
		1642	//.. elastic_predictor_stress.reportshort("ELASTIC PREDICTOR stress\n");
		1643	E_Pred_EPS.setStress( elastic_predictor_stress );
		1644
		1645	// double f_start = Criterion.f(start_stress);
		1646	// ::printf("SEMI BE############## f_start = %.10e\n",f_start);
		1647	double f_pred = getYS()->f( &E_Pred_EPS );
365	jeremic	1648	//::printf("SEMI BE############## f_pred = %.10e\n",f_pred);
13	jeremic	1649
		1650	// second of alternative predictors as seen in MAC page 170-171
		1651	if ( f_pred >= 0.0 )
		1652	{
		1653	//el_or_pl_range(1); // set to 1 ( plastic range )
		1654	// PREDICTOR FASE
		1655	//.. ::printf("\n\npredictor part step_counter = %d\n\n", step_counter);
		1656
		1657	dFods = getYS()->dFods( &E_Pred_EPS );
		1658	dQods = getPS()->dQods( &E_Pred_EPS );
		1659	//.... dFods.print("a","dF/ds ");
		1660	//.... dQods.print("a","dQ/ds ");
		1661
		1662	temp2 = (dFods("ij")E("ijkl"))dQods("kl");
		1663	temp2.null_indices();
		1664	lower = temp2.trace();
		1665
		1666	// Evaluating the hardening modulus: sum of (df/dq) qbar
		1667
		1668	//Of scalar internal var
		1669	hardMod_ = 0.0;
		1670
		1671	//Of 1st scalar internal vars
		1672	if ( getELS1() ) {
		1673	h_s[0] = getELS1()->h_s( &E_Pred_EPS, getPS());
		1674	xi_s[0] = getYS()->xi_s1( &E_Pred_EPS );
		1675	hardMod_ = hardMod_ + h_s[0] * xi_s[0];
		1676	}
		1677
		1678	//Of 2nd scalar internal vars
		1679	if ( getELS2() ) {
		1680	h_s[1] = getELS2()->h_s( &E_Pred_EPS, getPS());
		1681	xi_s[1] = getYS()->xi_s2( &E_Pred_EPS );
		1682	hardMod_ = hardMod_ + h_s[1] * xi_s[1];
		1683	}
		1684
		1685	//Of 3rd scalar internal vars
		1686	if ( getELS3() ) {
		1687	h_s[2] = getELS3()->h_s(&E_Pred_EPS, getPS());
		1688	xi_s[2] = getYS()->xi_s3( &E_Pred_EPS );
		1689	hardMod_ = hardMod_ + h_s[2] * xi_s[2];
		1690	}
		1691
		1692	//Of 4th scalar internal vars
		1693	if ( getELS4() ) {
		1694	h_s[3] = getELS4()->h_s( &E_Pred_EPS, getPS());
		1695	xi_s[3] = getYS()->xi_s4( &E_Pred_EPS );
		1696	hardMod_ = hardMod_ + h_s[3] * xi_s[3];
		1697	}
		1698
		1699	//Of tensorial internal var
		1700	// 1st tensorial var
		1701	if ( getELT1() ) {
		1702	h_t[0] = getELT1()->h_t(&E_Pred_EPS, getPS());
		1703	xi_t[0] = getYS()->xi_t1( &E_Pred_EPS );
		1704	tensor hm = (h_t[0])("ij") * (xi_t[0])("ij");
		1705	hardMod_ = hardMod_ + hm.trace();
		1706	}
		1707
		1708	// 2nd tensorial var
		1709	if ( getELT2() ) {
		1710	h_t[1] = getELT2()->h_t(&E_Pred_EPS, getPS());
		1711	xi_t[1] = getYS()->xi_t2( &E_Pred_EPS );
		1712	tensor hm = (h_t[1])("ij") * (xi_t[1])("ij");
		1713	hardMod_ = hardMod_ + hm.trace();
		1714	}
		1715
		1716	// 3rd tensorial var
		1717	if ( getELT3() ) {
		1718	h_t[2] = getELT3()->h_t(&E_Pred_EPS, getPS());
		1719	xi_t[2] = getYS()->xi_t3( &E_Pred_EPS );
		1720	tensor hm = (h_t[2])("ij") * (xi_t[2])("ij");
		1721	hardMod_ = hardMod_ + hm.trace();
		1722	}
		1723
		1724	// 4th tensorial var
		1725	if ( getELT4() ) {
		1726	h_t[3] = getELT4()->h_t(&E_Pred_EPS, getPS());
		1727	xi_t[3] = getYS()->xi_t4( &E_Pred_EPS );
		1728	tensor hm = (h_t[3])("ij") * (xi_t[3])("ij");
		1729	hardMod_ = hardMod_ + hm.trace();
		1730	}
		1731
		1732	// Subtract accumulated hardMod_ from lower
		1733	lower = lower - hardMod_;
		1734
		1735	//h_s = material_point.ELS1->h_s( &E_Pred_EPS, material_point.PS );
		1736	//xi_s = material_point.YS->xi_s1( &E_Pred_EPS );
		1737	//hardMod_ = h_s * xi_s;
		1738	//lower = lower - hardMod_;
		1739
		1740	////Of tensorial internal var
		1741	//h_t = material_point.ELT1->h_t(&E_Pred_EPS, material_point.PS);
		1742	//xi_t = material_point.YS->xi_t1( &E_Pred_EPS );
		1743	//tensor hm = h_t("ij") * xi_t("ij");
		1744	//hardMod_ = hm.trace();
		1745	//lower = lower - hardMod_;
		1746
		1747
		1748	Delta_lambda = f_pred/lower;
		1749	if ( Delta_lambda < 0.0 )
		1750	{
365	jeremic	1751	//::fprintf(stderr,"\nP!\n");
13	jeremic	1752	}
		1753	plastic_stress = (E("ijkl")dQods("kl"))(-Delta_lambda);
		1754	plastic_stress.null_indices();
		1755	//.. plastic_stress.reportshort("plastic stress predictor II\n");
		1756	//.. elastic_predictor_stress.reportshort("elastic predictor stress \n");
		1757	elastic_plastic_stress = elastic_predictor_stress + plastic_stress;
		1758	elastic_plastic_stress.null_indices();
		1759
		1760	SemibackwardEPS.setStress( elastic_plastic_stress );
		1761
		1762	////q_ast_entry = Criterion.kappa_get(elastic_plastic_stress);
		1763	//S = SemibackwardEPS.getScalarVar(1); // Get the old value of the internal var
		1764	//h_s = material_point.ELS1->h_s( &SemibackwardEPS, material_point.PS );
		1765	//dS = Delta_lambda * h_s ; // Increment to the internal scalar var
		1766	//h_t = material_point.ELT1->h_t( &SemibackwardEPS, material_point.PS );
		1767	//dT = Delta_lambda * h_t ; // Increment to the internal tensorial var
		1768
		1769	//Evolve the surfaces and hardening vars
		1770	int NS = SemibackwardEPS.getNScalarVar();
		1771	int NT = SemibackwardEPS.getNTensorVar();
		1772
458	jeremic	1773	for (int ii = 1; ii <= NS; ii++) {
13	jeremic	1774	dS = Delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var
		1775	S = SemibackwardEPS.getScalarVar(ii); // Get the old value of the scalar internal var
		1776	SemibackwardEPS.setScalarVar(ii, S + dS ); // Update internal scalar var
		1777	}
		1778
		1779
		1780	////current_lambda_set(Delta_lambda);
		1781	////q_ast = q_ast_entry + Dq_ast;
		1782	////Criterion.kappa_set( elastic_plastic_stress, q_ast);
		1783	//SemibackwardEPS.setScalarVar(1, S + dS );
		1784
		1785	//stresstensor new_T = T + dT;
		1786	//SemibackwardEPS.setTensorVar(1, new_T );
		1787
		1788	stresstensor new_T;
		1789
458	jeremic	1790	for (int ii = 1; ii <= NT; ii++) {
157	jeremic	1791	dT = h_t[ii-1] * Delta_lambda; // Increment to the tensor internal var
13	jeremic	1792	T = SemibackwardEPS.getTensorVar(ii); // Get the old value of the tensor internal var
		1793	new_T = T + dT;
		1794	SemibackwardEPS.setTensorVar(ii, new_T );
		1795	}
		1796
		1797
		1798	//return elastic_plastic_stress;
		1799	return SemibackwardEPS;
		1800	}
		1801	return E_Pred_EPS;
		1802	}
		1803
		1804
		1805
		1806
		1807	//================================================================================
		1808	// New EPState using Backward Euler Algorithm
		1809	//================================================================================
		1810	EPState Template3Dep::BackwardEulerEPState( const straintensor &strain_increment)
		1811
		1812	{
		1813	// Temp matertial point
		1814	//NDMaterial MP( material_point );
		1815	//NDMaterial *MP = this->getCopy();
130	jeremic	1816
		1817	// Volumetric strain
		1818	double st_vol = strain_increment.p_hydrostatic();
13	jeremic	1819
		1820	//EPState to be returned, it can be elastic or elastic-plastic EPState
		1821	EPState backwardEPS( * (this->getEPS()) );
		1822
		1823	EPState startEPS( *(this->getEPS()) );
		1824	stresstensor start_stress = startEPS.getStress();
		1825
		1826	//Output for plotting
		1827	cout.precision(5);
		1828	cout.width(10);
130	jeremic	1829	//cout << " strain_increment " << strain_increment << "\n";
13	jeremic	1830
		1831	cout.precision(5);
		1832	cout.width(10);
130	jeremic	1833	//cout << "start_stress " << start_stress;
13	jeremic	1834
		1835	// Pulling out some tensor and double definitions
		1836	tensor I2("I", 2, def_dim_2);
		1837	tensor I_ijkl("I", 4, def_dim_4);
		1838	I_ijkl = I2("ij")*I2("kl");
		1839	I_ijkl.null_indices();
		1840	tensor I_ikjl("I", 4, def_dim_4);
		1841	I_ikjl = I_ijkl.transpose0110();
		1842
		1843	//double Ey = Criterion.E();
		1844	//double nu = Criterion.nu();
		1845	//tensor E = StiffnessTensorE(Ey,nu);
		1846	tensor E = ElasticStiffnessTensor();
		1847	// tensor D = ComplianceTensorD(Ey,nu);
		1848	// stresstensor REAL_start_stress = start_stress;
		1849
		1850	tensor dFods(2, def_dim_2, 0.0);
		1851	tensor dQods(2, def_dim_2, 0.0);
		1852	// tensor dQodsextended(2, def_dim_2, 0.0);
		1853	// tensor d2Qodqast(2, def_dim_2, 0.0);
		1854	tensor temp2(2, def_dim_2, 0.0);
		1855	double lower = 0.0;
		1856	double Delta_lambda = 0.0; // Lambda
		1857	double delta_lambda = 0.0; // dLambda
		1858
		1859	double Felplpredictor = 0.0;
		1860	//Kai double absFelplpredictor = 0.0;
		1861	// double Ftolerance = pow(d_macheps(),(1.0/2.0))*1000000.00; //FORWARD no iterations
		1862	//double Ftolerance = pow( d_macheps(), 0.5)*1.00;
		1863
365	jeremic	1864	double Ftolerance = pow( d_macheps(), 0.5)*100000.00; //Zhaohui
13	jeremic	1865
		1866	// double Ftolerance = pow(d_macheps(),(1.0/2.0))*1.0;
		1867	// double entry_kappa_cone = Criterion.kappa_cone_get();
		1868	// double entry_kappa_cap = Criterion.kappa_cap_get();
		1869
		1870	tensor aC(2, def_dim_2, 0.0);
		1871	stresstensor BEstress;
		1872	stresstensor residual;
		1873	tensor d2Qoverds2( 4, def_dim_4, 0.0);
		1874	tensor T( 4, def_dim_4, 0.0);
		1875	tensor Tinv( 4, def_dim_4, 0.0);
		1876
		1877	double Fold = 0.0;
		1878	tensor temp3lower;
		1879	tensor temp5;
		1880	double temp6 = 0.0;
		1881	double upper = 0.0;
		1882
		1883	stresstensor dsigma;
		1884	//stresstensor Dsigma;
		1885	stresstensor sigmaBack;
		1886	straintensor dPlasticStrain; // delta plastic strain
		1887	straintensor PlasticStrain; // Total plastic strain
		1888
		1889	//double dq_ast = 0.0; // iterative change in internal variable (kappa in this case)
		1890	//double Dq_ast = 0.0; // incremental change in internal variable (kappa in this case)
		1891	//double q_ast = 0.0; // internal variable (kappa in this case)
		1892	//double q_ast_entry = 0.0; //internal variable from previous increment (kappa in this case)
		1893
		1894	int step_counter = 0;
		1895	//int MAX_STEP_COUNT = ;
		1896	// int MAX_STEP_COUNT = 0;
		1897	int flag = 0;
		1898
		1899	straintensor strain_incr = strain_increment;
		1900	strain_incr.null_indices();
		1901	stresstensor stress_increment = E("ijkl")*strain_incr("kl");
		1902	stress_increment.null_indices();
		1903
		1904	stresstensor Return_stress; // stress to be returned can be either predictor
		1905	// or elastic plastic stress.
		1906
		1907	EPState ElasticPredictorEPS( startEPS );
		1908	stresstensor elastic_predictor_stress = start_stress + stress_increment;
130	jeremic	1909
13	jeremic	1910	ElasticPredictorEPS.setStress( elastic_predictor_stress );
		1911	// elastic_predictor_stress.reportshortpqtheta("\n . . . . ELASTIC PREDICTOR stress");
		1912
130	jeremic	1913	cout.precision(5);
		1914	cout.width(10);
		1915	//cout << "elastic predictor " << elastic_predictor_stress << endln;
		1916
13	jeremic	1917	stresstensor elastic_plastic_predictor_stress;
		1918	EPState EP_PredictorEPS( startEPS );
		1919
		1920	//double f_start = material_point.YS->f( &startEPS );
		1921	//cout << " ************** f_start = " << f_start;
		1922	//::fprintf(stdout,"tst############## f_start = %.10e\n",f_start);
		1923	// f_pred = Criterion.f(elastic_predictor_stress);
		1924	//::fprintf(stdout,"tst############## f_pred = %.10e\n",f_pred);
		1925	double f_pred = getYS()->f( &ElasticPredictorEPS );
130	jeremic	1926	//cout << " ** f_pred ** " << f_pred << endln;
13	jeremic	1927	//int region = 5;
		1928
		1929	//double h_s = 0.0;
		1930	//double xi_s = 0.0;
		1931	double hardMod_ = 0.0;
		1932
		1933	double h_s[4] = {0.0, 0.0, 0.0, 0.0};
		1934	double xi_s[4] = {0.0, 0.0, 0.0, 0.0};
		1935	stresstensor h_t[4];
		1936	stresstensor xi_t[4];
		1937
		1938	//region
		1939	//check for the region than proceede
		1940	//region = Criterion.influence_region(elastic_predictor_stress);
		1941	//if ( region == 1 ) // apex gray region
		1942	// {
		1943	// double pc_ = pc();
		1944	// elastic_plastic_predictor_stress =
		1945	// elastic_plastic_predictor_stress.pqtheta2stress(pc_, 0.0, 0.0);
		1946	// return elastic_plastic_predictor_stress;
		1947	// }
		1948
		1949	if ( f_pred <= Ftolerance )
		1950	{
		1951
		1952	//Updating elastic strain increment
		1953	straintensor estrain = ElasticPredictorEPS.getElasticStrain();
		1954	straintensor tstrain = ElasticPredictorEPS.getStrain();
		1955	estrain = estrain + strain_increment;
		1956	tstrain = tstrain + strain_increment;
		1957	ElasticPredictorEPS.setElasticStrain( estrain );
		1958	ElasticPredictorEPS.setStrain( tstrain );
		1959	ElasticPredictorEPS.setdElasticStrain( strain_increment );
365	jeremic	1960	//cout<< "Elastic: Total strain" << tstrain << endl;
13	jeremic	1961
130	jeremic	1962	//if ( getELT1() )
		1963	// getELT1()->updateEeDm(&ElasticPredictorEPS, st_vol, 0.0);
		1964
13	jeremic	1965	//Set Elasto-Plastic stiffness tensor
		1966	ElasticPredictorEPS.setEep(E);
		1967	ElasticPredictorEPS.setConverged(TRUE);
130	jeremic	1968	//E.printshort(" BE -- Eep ");
13	jeremic	1969
		1970	backwardEPS = ElasticPredictorEPS;
		1971	return backwardEPS;
		1972	//cout << "\nbackwardEPS" << backwardEPS;
		1973	//cout << "\nElasticPredictorEPS " << ElasticPredictorEPS;
		1974
		1975	}
		1976	if ( f_pred > 0.0 )
		1977	{
		1978
		1979	//el_or_pl_range(1); // set to 1 ( plastic range )
		1980	//PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
		1981	//elastic_plastic_predictor_stress = Criterion.PredictorStress( start_stress,
		1982	// strain_increment,
		1983	// Criterion);
		1984
		1985	//Starting point by applying one Forward Euler step
		1986	EP_PredictorEPS = PredictorEPState( strain_incr);
		1987
		1988
		1989	//cout << " ----------Predictor Stress" << EP_PredictorEPS.getStress();
		1990	//Setting the starting EPState with the starting internal vars in EPState
		1991
		1992	//MP->setEPS( EP_PredictorEPS );
		1993
		1994	//PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
		1995	// IZBACIO ga 27 jan 97 . . .
		1996	// Delta_lambda = Criterion.current_lambda_get() ;
		1997	//::printf(" Delta_lambda = Criterion.current_lambda_get = %.8e\n", Delta_lambda);
		1998	//PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
		1999	//Felplpredictor = Criterion.f(elastic_plastic_predictor_stress);
		2000
		2001	Felplpredictor = getYS()->f(&EP_PredictorEPS);
		2002	//cout << " F_elplpredictor " << Felplpredictor << endln;
		2003
		2004
		2005	//Kai absFelplpredictor = fabs(Felplpredictor);
		2006	if ( fabs(Felplpredictor) <= Ftolerance )
		2007	{
		2008	backwardEPS = EP_PredictorEPS;
		2009
		2010	//Return_stress = elastic_plastic_predictor_stress;
		2011	flag = 1;
		2012	// this part should have been done up in "Criterion.PredictorStress" function.
		2013	// Criterion.kappa_set( Return_stress, q_ast) ;
		2014	// return Return_stress;
		2015	}
		2016	else {
		2017	////check for the region again and proceede
		2018	//// probably to see if predictor took you to any of the bad regions
		2019	//region = Criterion.influence_region(elastic_plastic_predictor_stress);
		2020	//if ( region == 1 ) // apex gray region
		2021	//{
		2022	// double pc_ = pc();
		2023	// elastic_plastic_predictor_stress =
		2024	// elastic_plastic_predictor_stress.pqtheta2stress(pc_, 0.0, 0.0);
		2025	// return elastic_plastic_predictor_stress;
		2026	//}
		2027
		2028	// NEWTON-RAPHSON RETURN
		2029	// residual stresstensor
		2030	//aC = Criterion.dQods(elastic_plastic_predictor_stress);
		2031	//dFods = Criterion.dFods(elastic_plastic_predictor_stress);
		2032	//dQods = Criterion.dQods(elastic_plastic_predictor_stress);
		2033	aC = getPS()->dQods( &EP_PredictorEPS );
		2034	dFods = getYS()->dFods( &EP_PredictorEPS );
		2035	dQods = getPS()->dQods( &EP_PredictorEPS );
		2036
		2037	// aC = Criterion.dQods(elastic_predictor_stress);
		2038	// dFods = Criterion.dFods(elastic_predictor_stress);
		2039	// dQods = Criterion.dQods(elastic_predictor_stress);
		2040
		2041
		2042	temp2 = (dFods("ij")E("ijkl"))dQods("kl");
		2043	temp2.null_indices();
		2044	lower = temp2.trace();
		2045
		2046	// Delta_lambda = f_pred/lower; //?????????
		2047	//::printf(" Delta_lambda = f_pred/lower = %.8e\n", Delta_lambda);
		2048	//// Delta_lambda = Felplpredictor/lower;
		2049	////::printf(" Delta_lambda = Felplpredictor/lower =%.8e \n", Delta_lambda);
		2050
		2051	// Original segment
130	jeremic	2052	//elastic_plastic_predictor_stress = elastic_predictor_stress - E("ijkl")aC("kl")Delta_lambda;
		2053	//EP_PredictorEPS.setStress( elastic_plastic_predictor_stress );
13	jeremic	2054
		2055	//Zhaohui modified, sometimes give much better convergence rate
130	jeremic	2056	elastic_plastic_predictor_stress = EP_PredictorEPS.getStress();
365	jeremic	2057	//cout << "elastic_plastic_predictor_stress" << elastic_plastic_predictor_stress;
13	jeremic	2058
		2059	cout.precision(5);
		2060	cout.width(10);
365	jeremic	2061	//cout << " " << EP_PredictorEPS.getStress().p_hydrostatic() << " ";
13	jeremic	2062
		2063	cout.precision(5);
		2064	cout.width(10);
365	jeremic	2065	//cout << EP_PredictorEPS.getStress().q_deviatoric()<< " ";
13	jeremic	2066
		2067	cout.precision(5);
		2068	cout.width(10);
365	jeremic	2069	//cout << Delta_lambda << endln;
13	jeremic	2070
		2071	//elastic_plastic_predictor_stress.reportshort("......elastic_plastic_predictor_stress");
		2072	//::printf(" F(elastic_plastic_predictor_stress) = %.8e\n", Criterion.f(elastic_plastic_predictor_stress));
		2073
		2074	//h_s = MP->ELS1->h_s( &EP_PredictorEPS, MP->PS );
		2075	////h_ = h(elastic_plastic_predictor_stress);
		2076	//// Dq_ast = Criterion.kappa_get(elastic_plastic_predictor_stress);
		2077
		2078	//q_ast_entry = Criterion.kappa_get(elastic_plastic_predictor_stress);
		2079	//Dq_ast = Delta_lambda * h_ * just_this_PP;
		2080	//q_ast = q_ast_entry + Dq_ast;
		2081
		2082	//Zhaohui comments: internal vars are alreary evolued in ForwardEulerEPS(...), not necessary here!!!
		2083	//..dS = Delta_lambda * h_s ; // Increment to the internal var
		2084	//..S = EP_PredictorEPS.getScalarVar(1); // Get the old value of the internal var
		2085	//..new_S = S + dS;
		2086	//..cout << "Internal Var : " << new_S << endln;
		2087	//..EP_PredictorEPS.setScalarVar(1, new_S); // Get the old value of the internal var
		2088
		2089	//::fprintf(stdout," ######## predictor --> Dq_ast=%.8e q_ast=%.8e\n", Dq_ast, q_ast);
		2090	//::fprintf(stderr," ######## predictor --> Dq_ast=%.8e q_ast=%.8e\n", Dq_ast, q_ast);
		2091
		2092	//Criterion.kappa_set( sigmaBack, q_ast); //????
		2093	//current_lambda_set(Delta_lambda); //????
		2094
		2095	//::printf(" Delta_lambda = %.8e\n", Delta_lambda);
		2096	//::printf("step = pre iteracija #############################-- q_ast = %.10e \n", q_ast);
		2097	//::printf("step = pre iteracija posle predictora ###########-- Dq_ast = %.10e \n",Dq_ast);
		2098	//**********
		2099	//**********
		2100	//::printf("\nDelta_lambda before BE = %.10e \n", Delta_lambda );
		2101	}
		2102
		2103	//========================== main part of iteration =======================
		2104	// while ( absFelplpredictor > Ftolerance &&
		2105
		2106	while ( fabs(Felplpredictor) > Ftolerance && step_counter < MAX_STEP_COUNT ) // if more iterations than prescribed
		2107	//out07may97 do
		2108	{
		2109	BEstress = elastic_predictor_stress - E("ijkl")aC("kl")Delta_lambda;
		2110	//BEstress.reportshort("......BEstress ");
		2111	///// BEstress = elastic_plastic_predictor_stress - E("ijkl")aC("kl")Delta_lambda;
		2112	BEstress.null_indices();
		2113	// Felplpredictor = Criterion.f(BEstress);
		2114	// ::printf("\nF_backward_Euler BE = %.10e \n", Felplpredictor);
		2115	residual = elastic_plastic_predictor_stress - BEstress;
		2116	//residual.reportshortpqtheta("\n......residual ");
		2117	// double ComplementaryEnergy = (residual("ij")D("ijkl")residual("ij")).trace();
		2118	//::printf("\n Residual ComplementaryEnergy = %.16e\n", ComplementaryEnergy);
		2119
		2120	///// residual = elastic_predictor_stress - BEstress;
		2121
		2122	//d2Qoverds2 = Criterion.d2Qods2(elastic_plastic_predictor_stress);
		2123	d2Qoverds2 = getPS()->d2Qods2( &EP_PredictorEPS );
		2124	//d2Qoverds2.print();
		2125
		2126	T = I_ikjl + E("ijkl")d2Qoverds2("klmn")Delta_lambda;
		2127	T.null_indices();
		2128
		2129	Tinv = T.inverse();
		2130	//dFods = Criterion.dFods(elastic_plastic_predictor_stress);
		2131	//dQods = Criterion.dQods(elastic_plastic_predictor_stress);
		2132	dFods = getYS()->dFods( &EP_PredictorEPS );
		2133	dQods = getPS()->dQods( &EP_PredictorEPS );
		2134
		2135	//Fold = Criterion.f(elastic_plastic_predictor_stress);
		2136	Fold = getYS()->f( &EP_PredictorEPS );
		2137
		2138	lower = 0.0; // this is old temp variable used here again :-)
		2139	//h_ = h(elastic_plastic_predictor_stress);
		2140	//xi_ = xi(elastic_plastic_predictor_stress);
		2141
		2142	//h_s = MP->ELS1->h_s( &EP_PredictorEPS, MP->PS );
		2143	//xi_s = MP->YS->xi_s1( &EP_PredictorEPS );
		2144	//hardMod_ = h_s * xi_s;
		2145
		2146	// Evaluating the hardening modulus: sum of (df/dq) qbar
		2147
		2148	hardMod_ = 0.0;
		2149	//Of 1st scalar internal vars
		2150	if ( getELS1() ) {
		2151	h_s[0] = getELS1()->h_s( &EP_PredictorEPS, getPS());
		2152	xi_s[0] = getYS()->xi_s1( &EP_PredictorEPS );
		2153	hardMod_ = hardMod_ + h_s[0] * xi_s[0];
		2154	}
		2155
		2156	//Of 2nd scalar internal vars
		2157	if ( getELS2() ) {
		2158	h_s[1] = getELS2()->h_s( &EP_PredictorEPS, getPS());
		2159	xi_s[1] = getYS()->xi_s2( &EP_PredictorEPS );
		2160	hardMod_ = hardMod_ + h_s[1] * xi_s[1];
		2161	}
		2162
		2163	//Of 3rd scalar internal vars
		2164	if ( getELS3() ) {
		2165	h_s[2] = getELS3()->h_s( &EP_PredictorEPS, getPS());
		2166	xi_s[2] = getYS()->xi_s3( &EP_PredictorEPS );
		2167	hardMod_ = hardMod_ + h_s[2] * xi_s[2];
		2168	}
		2169
		2170	//Of 4th scalar internal vars
		2171	if ( getELS4() ) {
		2172	h_s[3] = getELS4()->h_s( &EP_PredictorEPS, getPS());
		2173	xi_s[3] = getYS()->xi_s4( &EP_PredictorEPS );
		2174	hardMod_ = hardMod_ + h_s[3] * xi_s[3];
		2175	}
		2176
		2177	//Of tensorial internal var
		2178	// 1st tensorial var
		2179	if ( getELT1() ) {
		2180	h_t[0] = getELT1()->h_t( &EP_PredictorEPS, getPS());
		2181	xi_t[0] = getYS()->xi_t1( &EP_PredictorEPS );
		2182	tensor hm = (h_t[0])("ij") * (xi_t[0])("ij");
		2183	hardMod_ = hardMod_ + hm.trace();
		2184	}
		2185
		2186	// 2nd tensorial var
		2187	if ( getELT2() ) {
		2188	h_t[1] = getELT2()->h_t( &EP_PredictorEPS, getPS());
		2189	xi_t[1] = getYS()->xi_t2( &EP_PredictorEPS );
		2190	tensor hm = (h_t[1])("ij") * (xi_t[1])("ij");
		2191	hardMod_ = hardMod_ + hm.trace();
		2192	}
		2193
		2194	// 3rd tensorial var
		2195	if ( getELT3() ) {
		2196	h_t[2] = getELT3()->h_t( &EP_PredictorEPS, getPS());
		2197	xi_t[2] = getYS()->xi_t3( &EP_PredictorEPS );
		2198	tensor hm = (h_t[2])("ij") * (xi_t[2])("ij");
		2199	hardMod_ = hardMod_ + hm.trace();
		2200	}
		2201
		2202	// 4th tensorial var
		2203	if ( getELT4() ) {
		2204	h_t[3] = getELT4()->h_t( &EP_PredictorEPS, getPS());
		2205	xi_t[3] = getYS()->xi_t4( &EP_PredictorEPS );
		2206	tensor hm = (h_t[3])("ij") * (xi_t[3])("ij");
		2207	hardMod_ = hardMod_ + hm.trace();
		2208	}
		2209
		2210	// Subtract accumulated hardMod_ from lower
		2211	//lower = lower - hardMod_;
		2212
		2213	//hardMod_ = hardMod_ * just_this_PP;
		2214	//::printf("\n BackwardEulerStress .. hardMod_ = %.10e \n", hardMod_ );
		2215	//outfornow d2Qodqast = d2Qoverdqast(elastic_plastic_predictor_stress);
		2216	//outfornow dQodsextended = dQods + d2Qodqast * Delta_lambda * h_;
		2217	//outfornow temp3lower = dFods("mn")Tinv("ijmn")E("ijkl")*dQodsextended("kl");
458	jeremic	2218	temp3lower = dFods("mn")Tinv("ijmn")E("ijkl")*dQods("kl");
		2219	temp3lower.null_indices();
		2220	lower = temp3lower.trace();
13	jeremic	2221	lower = lower - hardMod_;
		2222
		2223	temp5 = (residual("ij") * Tinv("ijmn")) * dFods("mn");
		2224	temp6 = temp5.trace();
		2225	//The same as the above but more computation
		2226	//temp5 = dFods("mn") * residual("ij") * Tinv("ijmn");
		2227	//temp6 = temp5.trace();
		2228
		2229	upper = Fold - temp6;
		2230
		2231	//================================================================================
		2232	//dlambda
		2233	delta_lambda = upper / lower;
		2234	Delta_lambda = Delta_lambda + delta_lambda;
		2235
		2236	// Zhaohui_____10-01-2000 not sure xxxxxxxxxxxxxxxxxxx
		2237	dPlasticStrain = dQods("kl") * delta_lambda;
		2238	PlasticStrain = PlasticStrain + dPlasticStrain;
		2239
		2240	//::printf(" >> %d Delta_lambda = %.8e", step_counter, Delta_lambda);
		2241	// stari umesto dQodsextended za stari = dQods
		2242	//outfornow dsigma =
		2243	//outfornow ((residual("ij")*Tinv("ijmn"))+
		2244	//outfornow ((E("ijkl")dQodsextended("kl"))Tinv("ijmn")Delta_lambda) )-1.0;
		2245	//::printf(" Delta_lambda = %.8e\n", Delta_lambda);
458	jeremic	2246	dsigma = ( (residual("ij")*Tinv("ijmn") )+
		2247	( (E("ijkl")dQods("kl"))Tinv("ijmn")delta_lambda) )(-1.0); //*-1.0???
13	jeremic	2248
		2249	dsigma.null_indices();
		2250	//dsigma.reportshortpqtheta("\n......dsigma ");
		2251	//::printf(" ......... in NR loop delta_lambda = %.16e\n", delta_lambda);
		2252	//::printf(" ......... in NR loop Delta_lambda = %.16e\n", Delta_lambda);
		2253
		2254	//dq_ast = delta_lambda * h_ * just_this_PP;
		2255	//Dq_ast += dq_ast;
		2256
		2257	//dS = delta_lambda * h_s ; // Increment to the internal var
		2258	//S = EP_PredictorEPS.getScalarVar(1); // Get the old value of the internal var
		2259	//new_S = S + dS;
		2260	//EP_PredictorEPS.setScalarVar(1, new_S);
		2261
		2262	//Evolve the surfaces and hardening vars
		2263	int NS = EP_PredictorEPS.getNScalarVar();
		2264	int NT = EP_PredictorEPS.getNTensorVar();
		2265
		2266	double dS = 0;
		2267	double S = 0;
		2268	//double new_S = 0;
		2269
		2270	stresstensor dT;
		2271	stresstensor T;
		2272	stresstensor new_T;
458	jeremic	2273
		2274	for (int ii = 1; ii <= NS; ii++) {
13	jeremic	2275	dS = delta_lambda * h_s[ii-1] ; // Increment to the scalar internal var
		2276	S = EP_PredictorEPS.getScalarVar(ii); // Get the old value of the scalar internal var
		2277	EP_PredictorEPS.setScalarVar(ii, S + dS ); // Update internal scalar var
		2278	}
		2279
458	jeremic	2280	for (int ii = 1; ii <= NT; ii++) {
157	jeremic	2281	dT = h_t[ii-1] * delta_lambda; // Increment to the tensor internal var
13	jeremic	2282	T = EP_PredictorEPS.getTensorVar(ii); // Get the old value of the tensor internal var
		2283	new_T = T + dT;
		2284	EP_PredictorEPS.setTensorVar(ii, new_T ); // Update tensorial scalar var
130	jeremic	2285	}
13	jeremic	2286
		2287	//======= Dq_ast = Delta_lambda * h_ * just_this_PP;
		2288	//q_ast = q_ast_entry + Dq_ast;
		2289
		2290	//::fprintf(stdout," ######## step = %3d --> Dq_ast=%.8e q_ast=%.8e\n",
		2291	// step_counter, Dq_ast, q_ast);
		2292	//::fprintf(stderr," ######## step = %3d --> Dq_ast=%.8e q_ast=%.8e\n",
		2293	// step_counter, Dq_ast, q_ast);
		2294
		2295	//current_lambda_set(Delta_lambda);
		2296	//.... elastic_plastic_predictor_stress.reportshort("elplpredstress");
		2297	//.... dsigma.reportshort("dsigma");
		2298
		2299	//sigmaBack.reportshortpqtheta("\n before======== SigmaBack");
		2300	sigmaBack = elastic_plastic_predictor_stress + dsigma;
		2301	//sigmaBack.deviator().reportshort("\n after ======== SigmaBack");
		2302	//sigmaBack.reportshortpqtheta("\n after ======== SigmaBack");
		2303
130	jeremic	2304	//temp trick
		2305	if (sigmaBack.p_hydrostatic() > 0)
		2306	{
13	jeremic	2307	//====== sigmaBack = elastic_predictor_stress + Dsigma;
		2308	//sigmaBack.reportshortpqtheta("BE................ NR sigmaBack ");
		2309	//sigmaBack.reportAnim();
		2310	//::fprintf(stdout,"Anim BEpoint0%d = {Sin[theta]q, p, Cos[theta]q} \n",step_counter+1);
		2311	////::fprintf(stdout,"Anim BEpoint0%dP = Point[BEpoint0%d] \n",step_counter+1,step_counter+1);
		2312	//::fprintf(stdout,"Anim \n");
		2313
		2314	//Criterion.kappa_set( sigmaBack, q_ast) ;
		2315	EP_PredictorEPS.setStress( sigmaBack );
		2316
		2317	//Felplpredictor = Criterion.f(sigmaBack);
		2318	//Kai absFelplpredictor = fabs(Felplpredictor);
		2319	//::printf(" F_bE=%.10e (%.10e)\n", Felplpredictor,Ftolerance);
		2320	Felplpredictor = getYS()->f( &EP_PredictorEPS );
		2321	//::printf(" F_bE = %.10e (%.10e)\n", Felplpredictor, Ftolerance);
130	jeremic	2322	}
		2323	else
		2324	{
		2325	sigmaBack= sigmaBack.pqtheta2stress(0.1, 0.0, 0.0);
		2326	Felplpredictor = 0;
		2327	backwardEPS.setStress(sigmaBack);
		2328	backwardEPS.setConverged(TRUE);
		2329	return backwardEPS;
		2330	}
13	jeremic	2331
		2332	//double tempkappa1 = kappa_cone_get();
		2333	//double tempdFodeta = dFoverdeta(sigmaBack);
		2334	//double tempdetaodkappa = detaoverdkappa(tempkappa1);
		2335	//::printf(" h_=%.6e xi_=%.6e, dFodeta=%.6e, detaodkappa=%.6e, hardMod_=%.6e\n",
		2336	// h_, xi_,tempdFodeta, tempdetaodkappa, hardMod_);
		2337	//::printf(" upper = %.6e lower = %.6e\n", upper, lower);
		2338	//::printf(" q_ast_entry=%.6e Dq_ast=%.6e Delta_lambda=%.6e\n",
		2339	// q_ast_entry, Dq_ast, Delta_lambda);
		2340
		2341	// now prepare new step
		2342	elastic_plastic_predictor_stress = sigmaBack;
		2343
		2344	//Output for plotting
		2345	cout.precision(5);
		2346	cout.width(10);
365	jeremic	2347	//cout << " " << sigmaBack.p_hydrostatic() << " ";
130	jeremic	2348	//cerr << " " << sigmaBack.p_hydrostatic() << " ";
13	jeremic	2349
		2350	cout.precision(5);
		2351	cout.width(10);
365	jeremic	2352	//cout << sigmaBack.q_deviatoric() << " ";
130	jeremic	2353	//cerr << sigmaBack.q_deviatoric() << " ";
		2354
		2355	cout.precision(5);
		2356	cout.width(10);
365	jeremic	2357	//cout << " Felpl " << Felplpredictor;
130	jeremic	2358	//cerr << " Felpl " << Felplpredictor;
13	jeremic	2359
		2360	cout.precision(5);
		2361	cout.width(10);
365	jeremic	2362	//cout << " tol " << Ftolerance << endln;
130	jeremic	2363	//cerr << " tol " << Ftolerance << " " << step_counter << endln;
13	jeremic	2364
		2365	//::printf(" ........................... end of step %d\n", step_counter);// getchar();
		2366	step_counter++;
		2367	}
130	jeremic	2368	//cerr << " " << sigmaBack.p_hydrostatic() << " ";
		2369	//cerr << sigmaBack.q_deviatoric() << " ";
		2370	//cerr << " Felpl " << Felplpredictor;
		2371	//cerr << " tol " << Ftolerance << " " << step_counter << endln;
		2372
		2373	// Update E_Young and e according to current stress state before evaluate ElasticStiffnessTensor
		2374	int err = 0;
		2375	if ( getELT1() )
		2376	err = getELT1()->updateEeDm(&EP_PredictorEPS, st_vol, Delta_lambda);
		2377
		2378	//out07may97 while ( absFelplpredictor > Ftolerance &&
		2379	//out07may97 step_counter <= MAX_STEP_COUNT ); // if more iterations than prescribed
13	jeremic	2380
130	jeremic	2381	//**********
		2382	//**********
		2383	//**********
		2384	//**********
13	jeremic	2385	if ( step_counter >= MAX_STEP_COUNT )
		2386	{
130	jeremic	2387	//g3ErrorHandler->warning("Template3Dep::BackwardEulerEPState Step_counter > MAX_STEP_COUNT %d iterations", MAX_STEP_COUNT );
13	jeremic	2388	EP_PredictorEPS.setConverged( false );
		2389	//::exit(1);
		2390	}
		2391
		2392	// already set everything
		2393	if ( ( flag !=1) && (step_counter < MAX_STEP_COUNT) ) { // Need to genarate Eep and set strains and stresses
		2394
		2395	//Return_stress = elastic_plastic_predictor_stress;
		2396	//Criterion.kappa_set( Return_stress, q_ast) ;
		2397
		2398	// Generating Consistent Stiffness Tensor Eep
		2399	tensor I2("I", 2, def_dim_2);
		2400	tensor I_ijkl = I2("ij")*I2("kl");
		2401	I_ijkl.null_indices();
		2402	tensor I_ikjl = I_ijkl.transpose0110();
		2403
		2404
		2405	dQods = getPS()->dQods( &EP_PredictorEPS ); // this is m_ij
		2406	tensor temp2 = E("ijkl")*dQods("kl");
		2407	temp2.null_indices();
		2408	dFods = getYS()->dFods( &EP_PredictorEPS ); // this is n_ij
		2409	d2Qoverds2 = getPS()->d2Qods2( &EP_PredictorEPS );
		2410
		2411	tensor T = I_ikjl + E("ijkl")d2Qoverds2("klmn")Delta_lambda;
		2412	//tensor tt = E("ijkl")d2Qoverds2("klmn")Delta_lambda;
		2413	//tt.printshort("temp tt");
		2414	//T = I_ikjl + tt;
		2415	T.null_indices();
		2416	tensor Tinv = T.inverse();
		2417
		2418	tensor R = Tinv("ijmn")*E("ijkl");
		2419	R.null_indices();
		2420
		2421	// Evaluating the hardening modulus: sum of (df/dq) qbar at the final stress
		2422	hardMod_ = 0.0;
		2423	//Of 1st scalar internal vars
		2424	if ( getELS1() ) {
		2425	h_s[0] = getELS1()->h_s( &EP_PredictorEPS, getPS());
		2426	xi_s[0] = getYS()->xi_s1( &EP_PredictorEPS );
		2427	hardMod_ = hardMod_ + h_s[0] * xi_s[0];
		2428	}
		2429
		2430	//Of 2nd scalar internal vars
		2431	if ( getELS2() ) {
		2432	h_s[1] = getELS2()->h_s( &EP_PredictorEPS, getPS());
		2433	xi_s[1] = getYS()->xi_s2( &EP_PredictorEPS );
		2434	hardMod_ = hardMod_ + h_s[1] * xi_s[1];
		2435	}
		2436
		2437	//Of 3rd scalar internal vars
		2438	if ( getELS3() ) {
		2439	h_s[2] = getELS3()->h_s( &EP_PredictorEPS, getPS());
		2440	xi_s[2] = getYS()->xi_s3( &EP_PredictorEPS );
		2441	hardMod_ = hardMod_ + h_s[2] * xi_s[2];
		2442	}
		2443
		2444	//Of 4th scalar internal vars
		2445	if ( getELS4() ) {
		2446	h_s[3] = getELS4()->h_s( &EP_PredictorEPS, getPS());
		2447	xi_s[3] = getYS()->xi_s4( &EP_PredictorEPS );
		2448	hardMod_ = hardMod_ + h_s[3] * xi_s[3];
		2449	}
		2450
		2451	//Of tensorial internal var
		2452	// 1st tensorial var
		2453	if ( getELT1() ) {
		2454	h_t[0] = getELT1()->h_t( &EP_PredictorEPS, getPS());
		2455	xi_t[0] = getYS()->xi_t1( &EP_PredictorEPS );
		2456	tensor hm = (h_t[0])("ij") * (xi_t[0])("ij");
		2457	hardMod_ = hardMod_ + hm.trace();
		2458	}
		2459
		2460	// 2nd tensorial var
		2461	if ( getELT2() ) {
		2462	h_t[1] = getELT2()->h_t( &EP_PredictorEPS, getPS());
		2463	xi_t[1] = getYS()->xi_t2( &EP_PredictorEPS );
		2464	tensor hm = (h_t[1])("ij") * (xi_t[1])("ij");
		2465	hardMod_ = hardMod_ + hm.trace();
		2466	}
		2467
		2468	// 3rd tensorial var
		2469	if ( getELT3() ) {
		2470	h_t[2] = getELT3()->h_t( &EP_PredictorEPS, getPS());
		2471	xi_t[2] = getYS()->xi_t3( &EP_PredictorEPS );
		2472	tensor hm = (h_t[2])("ij") * (xi_t[2])("ij");
		2473	hardMod_ = hardMod_ + hm.trace();
		2474	}
		2475
		2476	// 4th tensorial var
		2477	if ( getELT4() ) {
		2478	h_t[3] = getELT4()->h_t( &EP_PredictorEPS, getPS());
		2479	xi_t[3] = getYS()->xi_t4( &EP_PredictorEPS );
		2480	tensor hm = (h_t[3])("ij") * (xi_t[3])("ij");
		2481	hardMod_ = hardMod_ + hm.trace();
		2482	}
		2483
		2484	tensor temp3lower = dFods("ot")R("otpq")dQods("pq");
		2485	temp3lower.null_indices();
		2486
		2487	double lower = temp3lower.trace();
		2488	lower = lower - hardMod_; // h
		2489
458	jeremic	2490	//tensor upper = R("pqkl")dQods("kl")dFods("ij")*R("ijmn");
		2491	tensor upper11 = R("pqkl")*dQods("kl");
		2492	tensor upper22 = dFods("ij")*R("ijmn");
		2493	upper11.null_indices();
		2494	upper22.null_indices();
		2495	tensor upper = upper11("pq")*upper22("mn");
		2496
		2497	upper.null_indices();
		2498	tensor Ep = upper*(1./lower);
13	jeremic	2499	tensor Eep = R - Ep; // elastoplastic constitutive tensor
		2500
		2501	//Set Elasto-Plastic stiffness tensor
		2502	EP_PredictorEPS.setEep(Eep);
		2503	EP_PredictorEPS.setConverged(TRUE);
		2504
		2505	//set plastic strain and total strain
		2506	straintensor elastic_strain = strain_increment - PlasticStrain; // elastic strain increment
		2507	straintensor estrain = EP_PredictorEPS.getElasticStrain(); //get old elastic strain
		2508	straintensor pstrain = EP_PredictorEPS.getPlasticStrain(); //get old plastic strain
		2509
		2510	straintensor tstrain = EP_PredictorEPS.getStrain(); //get old total strain
		2511	pstrain = pstrain + PlasticStrain;
		2512	estrain = estrain + elastic_strain;
		2513	tstrain = tstrain + strain_increment;
365	jeremic	2514	//cout<< "Plastic: Total strain" << tstrain <<endl;
13	jeremic	2515
		2516	//Setting de_p, de_e, total plastic, elastic strain, and total strain
		2517	EP_PredictorEPS.setdPlasticStrain( PlasticStrain );
		2518	EP_PredictorEPS.setdElasticStrain( elastic_strain );
		2519	EP_PredictorEPS.setPlasticStrain( pstrain );
		2520	EP_PredictorEPS.setElasticStrain( estrain );
		2521	EP_PredictorEPS.setStrain( tstrain );
		2522
		2523	backwardEPS = EP_PredictorEPS;
		2524	}
		2525
		2526	}
		2527	//return Return_stress;
		2528	backwardEPS = EP_PredictorEPS;
		2529	return backwardEPS;
		2530	}
		2531
		2532
		2533
		2534
		2535	//================================================================================
		2536	// New EPState using Forward Euler Subincrement Euler Algorithm
		2537	//================================================================================
130	jeremic	2538	EPState Template3Dep::FESubIncrementation( const straintensor & strain_increment,
13	jeremic	2539	int number_of_subincrements)
		2540	{
		2541
		2542	EPState old_EPS( *(this->getEPS()) );
		2543	EPState FESI_EPS( *(this->getEPS()) );
		2544	//NDMaterial MP( material_point );
		2545	//NDMaterial *MP = this->getCopy();
		2546	//cout << "in FESubIncrementation MP " << MP;
		2547
		2548	stresstensor back_stress;
		2549	stresstensor begin_stress = this->getEPS()->getStress();
		2550	//stresstensor begin_stress = start_stress;
		2551	//::fprintf(stderr,"{");
		2552
		2553	double sub = 1./( (double) number_of_subincrements );
		2554	//stresstensor elastic_subincremental_stress = stress_increment * sub;
		2555
130	jeremic	2556	straintensor tempp = strain_increment;
		2557	straintensor elastic_subincremental_strain = tempp * sub;
13	jeremic	2558	straintensor total_strain = elastic_subincremental_strain;
		2559	//elastic_subincremental_stress.reportshort("SUB INCREMENT in stresses\n");
365	jeremic	2560	//cout << "INCREMENT strain " << strain_increment << endln ;
		2561	//cout << "SUB INCREMENT strain " << elastic_subincremental_strain << endln ;
13	jeremic	2562
		2563	for( int steps=0 ; steps < number_of_subincrements ; steps++ ){
		2564
		2565	//start_stress.reportshort("START stress\n");
		2566	FESI_EPS = ForwardEulerEPState( elastic_subincremental_strain);
		2567
		2568	// Update the EPState in Template3Dep
		2569	this->setEPS( FESI_EPS );
		2570
		2571	back_stress = FESI_EPS.getStress();
		2572	cout.unsetf(ios::showpos);
		2573	cout << setw(4);
365	jeremic	2574	//cout << "Step No. " << steps << " ";
13	jeremic	2575
		2576	cout.setf(ios::scientific);
		2577	cout.setf(ios::showpos);
		2578	cout.precision(3);
		2579	cout << setw(7);
365	jeremic	2580	//cout << "p " << back_stress.p_hydrostatic() << " ";
		2581	//cout << setw(7);
		2582	//cout << "q " << back_stress.q_deviatoric() << " ";
		2583	//cout << setw(7);
		2584	//cout << " theta " << back_stress.theta() << " ";
		2585	//cout << setw(7);
		2586	//cout << "alfa1 " << FESI_EPS.getScalarVar(1) << " ";
		2587	//cout << setw(7);
		2588	//cout << "f = " << getYS()->f( &FESI_EPS ) << " "<< endln;
13	jeremic	2589
		2590	//begin_stress = back_stress;
		2591	//total_strain = total_strain + elastic_subincremental_strain;
		2592
		2593	}
		2594
		2595	// ::fprintf(stderr,"}");
		2596	this->setEPS( old_EPS );
		2597	return FESI_EPS;
		2598
		2599	}
		2600
		2601
		2602
		2603	//================================================================================
		2604	// New EPState using Backward Euler Subincrement Euler Algorithm
		2605	//================================================================================
		2606	EPState Template3Dep::BESubIncrementation( const straintensor & strain_increment,
		2607	int number_of_subincrements)
		2608	{
		2609	EPState old_EPS( *(this->getEPS()) );
		2610	EPState BESI_EPS( *(this->getEPS()) );
		2611	straintensor strain_incr = strain_increment;
		2612	//NDMaterial MP( material_point );
		2613	//NDMaterial *MP = getCopy();
		2614	//cout << "in FESubIncrementation MP " << MP;
		2615
		2616	stresstensor back_stress;
		2617	stresstensor begin_stress = getEPS()->getStress();
		2618	//stresstensor begin_stress = start_stress;
		2619	//::fprintf(stderr,"{");
		2620
		2621	double sub = 1./( (double) number_of_subincrements );
		2622	//stresstensor elastic_subincremental_stress = stress_increment * sub;
		2623
		2624	straintensor elastic_subincremental_strain = strain_incr * sub;
		2625	straintensor total_strain = elastic_subincremental_strain;
		2626	//elastic_subincremental_stress.reportshort("SUB INCREMENT in stresses\n");
		2627
		2628	for( int steps=0 ; steps < number_of_subincrements ; steps++ ){
		2629
		2630	//start_stress.reportshort("START stress\n");
		2631	BESI_EPS = BackwardEulerEPState( elastic_subincremental_strain);
130	jeremic	2632	//if ( !BESI_EPS.getConverged() )
13	jeremic	2633	this->setEPS( BESI_EPS );
130	jeremic	2634	//else {
		2635	//g3ErrorHandler->warning("Template3Dep::BESubIncrementation failed to converge using %d step sub-BackwardEuler Algor.", number_of_subincrements );
		2636	//exit(1);
		2637	//g3ErrorHandler->fatal("Template3Dep::BESubIncrementation failed to converge using %d step sub-BackwardEuler Algor.", number_of_subincrements );
		2638	//exit(1);
		2639	//}
13	jeremic	2640
		2641	//back_stress = BESI_EPS.getStress();
		2642	//cout.unsetf(ios::showpos);
		2643	//cout << setw(4);
		2644	//cout << "Step No. " << steps << " ";
		2645	//
		2646	//cout.setf(ios::scientific);
		2647	//cout.setf(ios::showpos);
		2648	//cout.precision(3);
		2649	//cout << setw(7);
		2650	//cout << "p " << back_stress.p_hydrostatic() << " ";
		2651	//cout << setw(7);
		2652	//cout << "q " << back_stress.q_deviatoric() << " ";
		2653	//cout << setw(7);
		2654	//cout << "alfa1 " << BESI_EPS.getScalarVar(1) << " ";
		2655	//cout << setw(7);
		2656	//cout << "f = " << MP->YS->f( &BESI_EPS ) << " "<< endln;
		2657
		2658	////begin_stress = back_stress;
		2659	////total_strain = total_strain + elastic_subincremental_strain;
		2660
		2661	}
		2662
		2663	// ::fprintf(stderr,"}");
		2664	this->setEPS( old_EPS );
		2665	return BESI_EPS;
		2666
		2667	}
		2668
		2669
		2670
		2671	////================================================================================
		2672	//// Routine used to generate elastic stiffness tensor E
		2673	////================================================================================
		2674	//tensor Template3Dep::ElasticStiffnessTensor( double E, double nu) const
		2675	// {
		2676	// tensor ret( 4, def_dim_4, 0.0 );
		2677	//
		2678	// tensor I2("I", 2, def_dim_2);
		2679	// tensor I_ijkl = I2("ij")*I2("kl");
		2680	// I_ijkl.null_indices();
		2681	// tensor I_ikjl = I_ijkl.transpose0110();
		2682	// tensor I_iljk = I_ijkl.transpose0111();
		2683	// tensor I4s = (I_ikjl+I_iljk)*0.5;
		2684	//
		2685	// // Building elasticity tensor
		2686	// ret = (I_ijkl((Enu2.0)/(2.0(1.0+nu)(1-2.0nu)))) + (I4s*(E/((1.0+nu))));
		2687	//
		2688	// return ret;
		2689	// }
		2690	//
		2691	//
		2692	////================================================================================
		2693	//// Routine used to generate elastic compliance tensor D
		2694	////================================================================================
		2695	//
		2696	//tensor Template3Dep::ElasticComplianceTensor( double E, double nu) const
		2697	// {
		2698	// if (E == 0) {
		2699	// cout << endln << "Ey = 0! Can't give you D!!" << endln;
		2700	// exit(1);
		2701	// }
		2702	//
		2703	// tensor ret( 4, def_dim_4, 0.0 );
		2704	// //tensor ret;
		2705	//
		2706	// tensor I2("I", 2, def_dim_2);
		2707	// tensor I_ijkl = I2("ij")*I2("kl");
		2708	// I_ijkl.null_indices();
		2709	// tensor I_ikjl = I_ijkl.transpose0110();
		2710	// tensor I_iljk = I_ijkl.transpose0111();
		2711	// tensor I4s = (I_ikjl+I_iljk)*0.5;
		2712	//
		2713	// // Building compliance tensor
		2714	// ret = (I_ijkl * (-nu/E)) + (I4s * ((1.0+nu)/E));
		2715	//
		2716	// return ret;
		2717	// }
		2718	//
		2719
		2720	//================================================================================
		2721	// trying to find intersection point
		2722	// according to M. Crisfield's book
		2723	// "Non-linear Finite Element Analysis of Solids and Structures "
		2724	// chapter 6.6.1 page 168.
		2725	//================================================================================
		2726	stresstensor Template3Dep::yield_surface_cross(const stresstensor & start_stress,
		2727	const stresstensor & end_stress)
		2728	{
		2729	// Bounds
		2730	double x1 = 0.0;
		2731	double x2 = 1.0;
		2732
		2733	// accuracy
		2734	double const TOL = 1.0e-9;
		2735	//cout << "start_stress "<< start_stress;
		2736	//cout << "end_stress " << end_stress;
		2737	//end_stress.reportshortpqtheta("end stress");
		2738
		2739	double a = zbrentstress( start_stress, end_stress, x1, x2, TOL ); // Defined below
		2740	// ::printf("\n**\n a = %lf \n**\n",a);
		2741
		2742	stresstensor delta_stress = end_stress - start_stress;
		2743	stresstensor intersection_stress = start_stress + delta_stress * a;
		2744	//*** intersection_stress.reportshort("FINAL Intersection stress\n");
		2745
		2746	return intersection_stress;
		2747
		2748	}
		2749
		2750
		2751	//================================================================================
		2752	// Routine used by yield_surface_cross to
		2753	// find the stresstensor at cross point
		2754	//================================================================================
		2755
		2756	double Template3Dep::zbrentstress(const stresstensor & start_stress,
		2757	const stresstensor & end_stress,
		2758	double x1, double x2, double tol)
		2759	{
		2760	double EPS = d_macheps();
		2761
		2762	int iter;
		2763	double a=x1;
		2764	double b=x2;
		2765	double c=0.0;
		2766	double d=0.0;
		2767	double e=0.0;
		2768	double min1=0.0;
		2769	double min2=0.0;
		2770	double fa=func(start_stress, end_stress, a);
		2771	double fb=func(start_stress, end_stress, b);
		2772	//cout << "fb = " << fb;
		2773
		2774	double fc=0.0;
		2775	double p=0.0;
		2776	double q=0.0;
		2777	double r=0.0;
		2778	double s=0.0;
		2779	double tol1=0.0;
		2780	double xm=0.0;
		2781	//::printf("\n############# zbrentstress iterations --\n");
		2782	if (fb*fa > 0.0)
		2783	{
		2784	::printf("\a\nRoot must be bracketed in ZBRENTstress\n");
		2785	::exit(1);
		2786	}
		2787	fc=fb;
		2788	for ( iter=1 ; iter<=ITMAX ; iter++ )
		2789	{
		2790	//::printf("iter No. = %d ; b = %16.10lf\n", iter, b);
		2791	if (fb*fc > 0.0)
		2792	{
		2793	c=a;
		2794	fc=fa;
		2795	e=d=b-a;
		2796	}
		2797	if (fabs(fc) < fabs(fb))
		2798	{
		2799	a=b;
		2800	b=c;
		2801	c=a;
		2802	fa=fb;
		2803	fb=fc;
		2804	fc=fa;
		2805	}
		2806	tol1=2.0EPSfabs(b)+0.5*tol;
		2807	xm=0.5*(c-b);
		2808	if (fabs(xm) <= tol1 \|\| fb == 0.0) return b;
		2809	if (fabs(e) >= tol1 && fabs(fa) > fabs(fb))
		2810	{
		2811	s=fb/fa;
		2812	if (a == c)
		2813	{
		2814	p=2.0xms;
		2815	q=1.0-s;
		2816	}
		2817	else
		2818	{
		2819	q=fa/fc;
		2820	r=fb/fc;
		2821	p=s(2.0xmq(q-r)-(b-a)*(r-1.0));
		2822	q=(q-1.0)(r-1.0)(s-1.0);
		2823	}
		2824	if (p > 0.0) q = -q;
		2825	p=fabs(p);
		2826	min1=3.0xmq-fabs(tol1*q);
		2827	min2=fabs(e*q);
		2828	if (2.0*p < (min1 < min2 ? min1 : min2))
		2829	{
		2830	e=d;
		2831	d=p/q;
		2832	}
		2833	else
		2834	{
		2835	d=xm;
		2836	e=d;
		2837	}
		2838	}
		2839	else
		2840	{
		2841	d=xm;
		2842	e=d;
		2843	}
		2844	a=b;
		2845	fa=fb;
		2846	if (fabs(d) > tol1)
		2847	b += d;
		2848	else
		2849	b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
		2850	fb=func(start_stress, end_stress, b);
		2851	}
365	jeremic	2852	//::printf("\a\nMaximum number of iterations exceeded in zbrentstress\n");
13	jeremic	2853	return 0.0; // this just to full the compiler because of the warnings
		2854	}
		2855
		2856
		2857	//================================================================================
		2858	// routine used by zbrentstress, takes an alfa and returns the
		2859	// yield function value for that alfa
		2860	//================================================================================
		2861	double Template3Dep::func(const stresstensor & start_stress,
		2862	const stresstensor & end_stress,
		2863	double alfa )
		2864	{
		2865
		2866	EPState *tempEPS = getEPS()->newObj();
		2867	stresstensor delta_stress = end_stress - start_stress;
		2868	stresstensor intersection_stress = start_stress + delta_stress*alfa;
		2869
		2870	tempEPS->setStress(intersection_stress);
		2871
		2872	//cout << "tempEPS" << tempEPS;
		2873
		2874	double f = getYS()->f( tempEPS );
		2875	return f;
		2876	}
		2877
		2878	//================================================================================
12	jeremic	2879	// Overloaded Insertion Operator
		2880	// prints an EPState's contents
		2881	//================================================================================
		2882	ostream& operator<< (ostream& os, const Template3Dep & MP)
		2883	{
		2884	os << endln << "Template3Dep: " << endln;
		2885	os << "\ttag: " << MP.getTag() << endln;
		2886	os << "=================================================================" << endln;
		2887	MP.getYS()->print();
		2888	MP.getPS()->print();
		2889	MP.getEPS()->print();
		2890
		2891	cout << endln << "Scalar Evolution Laws: " << endln;
		2892	if ( MP.ELS1 ){
130	jeremic	2893	os << "\nFor 1st scalar var:\n";
12	jeremic	2894	MP.ELS1->print();
		2895	}
		2896
		2897	if ( MP.ELS2 ){
130	jeremic	2898	os << "\nFor 2nd scalar var:\n";
12	jeremic	2899	MP.ELS2->print();
		2900	}
		2901
		2902	if ( MP.ELS3 ){
130	jeremic	2903	os << "\nFor 3rd scalar var:\n";
12	jeremic	2904	MP.ELS3->print();
		2905	}
		2906
		2907	if ( MP.ELS4 ){
130	jeremic	2908	os << "\nFor 4th scalar var:\n";
12	jeremic	2909	MP.ELS4->print();
		2910	}
		2911
		2912
		2913	cout << endln << "Tensorial Evolution Laws: " << endln;
		2914	if ( MP.ELT1 ){
130	jeremic	2915	os << "\nFor 1st tensorial var:\n";
12	jeremic	2916	MP.ELT1->print();
		2917	}
		2918	if ( MP.ELT2 ){
130	jeremic	2919	os << "\nFor 2nd tensorial var:\n";
12	jeremic	2920	MP.ELT2->print();
		2921	}
		2922	if ( MP.ELT3 ){
130	jeremic	2923	os << "\nFor 3rd tensorial var:\n";
12	jeremic	2924	MP.ELT3->print();
		2925	}
		2926	if ( MP.ELT4 ){
130	jeremic	2927	os << "\nFor 4th tensorial var:\n";
12	jeremic	2928	MP.ELT4->print();
		2929	}
		2930
		2931	os << endln;
		2932	return os;
		2933	}
		2934
458	jeremic	2935
12	jeremic	2936	#endif
		2937

Subversion Repositories OpenSees

(root)/trunk/SRC/material/nD/Template3Dep/Template3Dep.cpp @ 458 - Rev 458