WebSVN - OpenSees - Blame - Rev 4147 - /trunk/SRC/material/uniaxial/snap/Bilin.cpp

Rev	Author	Line No.	Line
4147	fmk	1	/* ****************************************************************
		2	OpenSees - Open System for Earthquake Engineering Simulation
		3	Pacific Earthquake Engineering Research Center
		4
		5
		6	(C) Copyright 1999, The Regents of the University of California
		7	All Rights Reserved.
		8
		9	Commercial use of this program without express permission of the
		10	University of California, Berkeley, is strictly prohibited. See
		11	file 'COPYRIGHT' in main directory for information on usage and
		12	redistribution, and for a DISCLAIMER OF ALL WARRANTIES.
		13
		14	Developed by:
		15	Frank McKenna (fmckenna@ce.berkeley.edu)
		16	Gregory L. Fenves (fenves@ce.berkeley.edu)
		17	Filip C. Filippou (filippou@ce.berkeley.edu)
		18
		19	**************************************************************** */
		20
		21
		22	#include <math.h>
		23
		24	#include <Bilin.h>
		25	#include <Vector.h>
		26	#include <Channel.h>
		27
		28	#include <OPS_Globals.h>
		29
		30	#define OPS_Export
		31
		32	static int numSAWSMaterials = 0;
		33
		34	OPS_Export void *
		35	OPS_NewSAWSMaterial()
		36	{
		37	if (numSAWSMaterials == 0) {
		38	numSAWSMaterials++;
		39	OPS_Error("Bilin unaxial material - Written by D. Lignos, Stanfurd 2009\n", 1);
		40	}
		41
		42	// Pointer to a uniaxial material that will be returned
		43	UniaxialMaterial *theMaterial = 0;
		44
		45	int iData[1];
		46	double dData[23];
		47	int numData = 1;
		48
		49	if (OPS_GetIntInput(&numData, iData) != 0) {
		50	opserr << "WARNING invalid uniaxialMaterial SAWSMaterial tag" << endln;
		51	return 0;
		52	}
		53
		54	numData = 23;
		55	if (OPS_GetDoubleInput(&numData, dData) != 0) {
		56	opserr << "Invalid Args want: uniaxialMaterial Bilin tag? Ke? As? AsNeg? My_pos? My_neg? LamdaS? ";
		57	opserr << "LamdaK? LamdaA? LamdaD? Cs? Ck? Ca? Cd? Thetap_pos? Thetap_neg? Thetapc_pos? Thetapc_neg?K? ";
		58	opserr << "KNeg? Thetau_pos? Thetau_neg? PDPlus? PDNeg\n";
		59	return 0;
		60	}
		61
		62	// Parsing was successful, allocate the material
		63	theMaterial = new Bilin(iData[0],
		64	dData[0], dData[1], dData[2], dData[3], dData[4],
		65	dData[5], dData[6], dData[7], dData[8], dData[9],
		66	dData[10], dData[11], dData[12], dData[13], dData[14],
		67	dData[15], dData[16], dData[17], dData[18], dData[19],
		68	dData[22], dData[21], dData[22]);
		69
		70	if (theMaterial == 0) {
		71	opserr << "WARNING could not create uniaxialMaterial of type Bilin Material\n";
		72	return 0;
		73	}
		74
		75	return theMaterial;
		76	}
		77
		78	Bilin::Bilin(int tag, double p_Ke,double p_As,double p_AsNeg,double p_My_pos,double p_My_neg,double p_LamdaS,
		79	double p_LamdaK,double p_LamdaA,double p_LamdaD,double p_Cs,double p_Ck,double p_Ca,double p_Cd,
		80	double p_Thetap_pos,double p_Thetap_neg,double p_Thetapc_pos,double p_Thetapc_neg,double p_K,double p_KNeg,
		81	double p_Thetau_pos,double p_Thetau_neg,double p_PDPlus,double p_PDNeg)
		82	:UniaxialMaterial(tag,MAT_TAG_Bilin),Ke(p_Ke), As(p_As), AsNeg(p_AsNeg), My_pos(p_My_pos), My_neg(p_My_neg),
		83	LamdaS(p_LamdaS), LamdaK(p_LamdaK),LamdaA(p_LamdaA),LamdaD(p_LamdaD), Cs(p_Cs), Ck(p_Ck), Ca(p_Ca),Cd(p_Cd),
		84	Thetap_pos(p_Thetap_pos), Thetap_neg(p_Thetap_neg), Thetapc_pos(p_Thetapc_pos),Thetapc_neg(p_Thetapc_neg),
		85	K(p_K), KNeg(p_KNeg),Thetau_pos(p_Thetau_pos), Thetau_neg(p_Thetau_neg), PDPlus(p_PDPlus), PDNeg(p_PDNeg)
		86	{
		87	//initialize variables
		88	this->revertToStart();
		89	//this->revertToLastCommit();
		90	}
		91
		92	Bilin::Bilin()
		93	:UniaxialMaterial(0,MAT_TAG_Bilin),
		94	Ke(0), As(0), AsNeg(0), My_pos(0), My_neg(0),
		95	LamdaS(0), LamdaK(0),LamdaA(0),LamdaD(0), Cs(0), Ck(0), Ca(0),Cd(0),
		96	Thetap_pos(0), Thetap_neg(0), Thetapc_pos(0),Thetapc_neg(0),
		97	K(0), KNeg(0),Thetau_pos(0), Thetau_neg(0), PDPlus(0), PDNeg(0)
		98	{
		99	this->revertToStart();
		100	}
		101
		102	Bilin::~Bilin()
		103	{
		104	// does nothing
		105	}
		106
		107	int
		108	Bilin::setTrialStrain(double strain, double strainRate)
		109	{
		110	//all variables to the last commit
		111	this->revertToLastCommit();
		112
		113	//Here I declare function variables
		114	double Ue,ddise,deltaD,d,temp_1_farzin,temp,betas,betak,betad,ekhard,
		115	dBoundPos,dBoundNeg,f1,f2,fn,e1,e2,a2,fNewLoadPos,
		116	f,xDevPos1,yDevPos1,xDevPos2,yDevPos2,xDevPos,yDevPos,fNewLoadNeg,xDevNeg1,
		117	yDevNeg1,xDevNeg2,yDevNeg2,xDevNeg,yDevNeg,Enrgi,v1,d1,ener;
		118
		119	int NoCollapse,kst,jcode;
		120
		121
		122	//state determination algorithm: defines the current force and tangent stiffness
		123	U=strain; //set trial displacement
		124
		125	//******************************matlab code***********************************************************
		126	//// total displacement
		127	Ue = U;
		128	//// incremental displacement
		129	ddise=Ue-CU; //the incremental displacement of the element (from the last commited)
		130	//// update Uprev
		131	Uprev=Ue;
		132	////
		133	//// calculation of deltaD
		134	deltaD=ddise;
		135	d=Ue;
		136	dP = d-deltaD;
		137	if (d>0.0)
		138	{
		139
		140	//function call
		141	interPoint(temp_1_farzin,temp,0.0,fCapRefPos,capSlopeelstk,0.0,ResfacfyieldPos,0.0);
		142	if (d<temp_1_farzin){
		143	iNoFpos = 0;
		144	LP=0;
		145	}
		146	} else {
		147	interPoint(temp_1_farzin,temp,0.0,fCapRefNeg,capSlopeNegelstk,0.0,ResfacNegfyieldNeg,0.0);
		148	if (d>temp_1_farzin)
		149	{
		150	iNoFneg = 0;
		151	LN=0;
		152	}
		153	}
		154	//// Other variables
		155	flagdeg = 0;
		156	betas = 0.0;
		157	betak = 0.0;
		158	betad = 0.0;
		159	//// yield siplacements
		160	dyieldPos = fyieldPos/elstk;
		161	dyieldNeg = fyieldNeg/elstk;
		162	//// Initialize parameters in the first cycle
		163	if (kon==0) {
		164	//Because I included the statement revertToLastCommit above:
		165	elstk = Ke;
		166	fyieldPos = My_pos;
		167	fyieldNeg = My_neg;
		168	alpha = As;
		169	alphaN = AsNeg;
		170	ecaps = LamdaS/(fyieldPos/(elstk));
		171	ecapk = LamdaK/(fyieldPos/(elstk));
		172	ecapa = LamdaA/(fyieldPos/(elstk));
		173	ecapd = LamdaD/(fyieldPos/(elstk));
		174	cs = Cs;
		175	ck = Ck;
		176	ca = Ca;
		177	cd = Cd;
		178	capDispPos = Thetap_pos+fyieldPos/elstk;
		179	capDispNeg = -Thetap_neg+fyieldNeg/elstk;
		180	Resfac = K;
		181	ResfacNeg = KNeg;
		182	myFcapping =-(fyieldPos+alphaelstkcapDispPos);
		183	capSlope = myFcapping/(Thetapc_pos*elstk);
		184	myFcappingNeg = -(abs(fyieldNeg)+alphaelstkabs(capDispNeg));
		185	capSlopeNeg = myFcappingNeg/(Thetapc_neg*elstk);
		186	fracDispPos = Thetau_pos;
		187	fracDispNeg = -Thetau_neg;
		188	DPlus = PDPlus;
		189	DNeg = PDNeg;
		190	stif = elstk;
		191	ekP = elstk;
		192	flagControlResponse = 0;
		193	Tangent=Ke;
		194	//stop
		195	ekhard = elstk*alpha;
		196	alphaNeg=alphaN;
		197	alphaPos=alpha;
		198	ekhardPos = ekhard;
		199	ekhardNeg = ekhard;
		200	dyieldPos = fyieldPos/elstk;
		201	dyieldNeg = fyieldNeg/elstk;
		202	Enrgts = fyieldPosdyieldPosecaps;
		203	Enrgtk = fyieldPosdyieldPosecapk;
		204	Enrgtd = fyieldPosdyieldPosecapd;
		205	dmax = dyieldPos;
		206	dmin = dyieldNeg;
		207	ekP = elstk;
		208	ekunload = elstk;
		209	ekexcurs = elstk;
		210	Enrgtot = 0.0;
		211	Enrgc = 0.0;
		212	fyPos = fyieldPos;
		213	fyNeg = fyieldNeg;
		214	dyNeg=dyieldNeg;
		215	dyPos=dyieldPos;
		216	resSn = fyieldPos;
		217	resSp = fyieldNeg;
		218	cpPos = capDispPos;
		219	cpNeg = capDispNeg;
		220	fPeakPos=fyieldPos+ekhard*(capDispPos-dyieldPos);
		221	fPeakNeg=fyieldNeg+ekhard*(capDispNeg-dyieldNeg);
		222	flagControlResponse = flagControlResponse; // This flag Controls the ultimate rotation
		223	DPlus = DPlus; // Slab Effect Positive Direction
		224	DNeg = DNeg; // Slab Effect Negative Direction
		225	if (cpPos<dyieldPos) {
		226	fPeakPos =fyieldPos*cpPos/dyieldPos;
		227	}
		228	if (cpNeg>dyieldNeg) {
		229	fPeakNeg =fyieldNeg*cpNeg/dyieldNeg;
		230	}
		231	fCapPos = fPeakPos;
		232	fCapNeg = fPeakNeg;
		233	LP = 0;
		234	LN = 0;
		235	fLimPos = 0;
		236	fLimNeg = 0;
		237	dLimPos = 0;
		238	dLimNeg = 0;
		239	dBoundPos = 0;
		240	dBoundNeg = 0;
		241	iNoFpos = 0;
		242	iNoFneg = 0;
		243	interup=0;
		244	fCapRefPos=-capSlopeelstkcapDispPos+fPeakPos;
		245	fCapRefNeg=-capSlopeNegelstkcapDispNeg+fPeakNeg;
		246	capSlopeOrig = capSlope;
		247	capSlopeOrigNeg = capSlopeNeg;
		248	flagstopdeg = 0;
		249	dCap2Pos = 0.0;
		250	dCap1Pos = 0.0;
		251	dCap2Neg = 0.0;
		252	dCap1Neg = 0.0;
		253	dyPos = 0.0;
		254	dyNeg = 0.0;
		255	f1 = 0.0;
		256	f2 = 0.0;
		257	fmin = 0.0;
		258	fmax = 0.0;
		259	RSE = 0.0;
		260	fn = 0.0;
		261	dn = 0.0;
		262	e1 = 0.0;
		263	e2 = 0.0;
		264	Enrgi = 0.0;
		265	ekt = 0.0;
		266	a2 = 0.0;
		267	NoCollapse = 0;
		268	iDeg = 0;
		269	if(deltaD>=0.0){
		270	kon = 1;
		271	} else {
		272	kon = 2;
		273	}
		274	kst = 1;
		275	}
		276	////
		277	// ***************** S T A R T S B I G L O O P **************
		278	// IF D E L T A > 0 - - - - - - - - - - - - - - - - - - - - - - - -
		279	if (deltaD>=0.0) {
		280	if (iNoFpos==1) {
		281	interPoint(dNewLoadPos,fNewLoadPos,dyNeg,fyNeg,ekhardNeg,0.0,0.0,0.0);
		282	}
		283	//If there is a corner changing delta from negative to positive
		284	if (kon==2) {
		285	kon = 1;
		286	dlstNeg = dP;
		287	flstNeg = fP;
		288
		289	if (dlstNeg<=dmin){
		290	fLimNeg = flstNeg;
		291	dLimNeg = dlstNeg;
		292	}
		293
		294	if (resSn>0){
		295	RSE = 0.5fPfP/ekunload;
		296	} else {
		297	RSE=0.5(fP+resSn)(sn-dP)+0.5resSnresSn/(elstk*alphaPos);
		298	}
		299
		300	a2 = Enrgtk-(Enrgtot-RSE);
		301	if((a2<=0.0) && (Enrgtk!=0.0)) {
		302	flgstop=1;
		303	}
		304
		305	if((ecapk!=0.0) && (d<sp) && abs(capSlope)>=1.0e-3 && (abs(capSlopeNeg)>=1.0e-3)){
		306	betak = pow(((Enrgc-RSE)/(Enrgtk-(Enrgtot-RSE))),ck);
		307	if(((Enrgtot-RSE)>=Enrgtk)\|\|(betak>=1.0)) {
		308	betak = 1.0;
		309	}
		310	if( flagstopdeg !=1 && flagdeg !=1) {
		311	ekunload = ekexcurs*(1-betak);
		312	} else {
		313	ekunload = ekexcurs;
		314	}
		315	if(ekunload<=0.1*elstk) {
		316	ekunload = 0.1*elstk;
		317	}
		318	}
		319
		320
		321	//// sn CALCULATION-------------------------------------------------
		322
		323	if((dmin<dyieldNeg)\|\|(dmax>dyieldPos)) {
		324	if((dP<sp)\|\|(((dP>sp)&&(ekP==ekunload)))) {
		325
		326	// call snCalc(sn,resSn,snEnv,resSnEnv,dP,fP,ekunload,alphaPos,dyPos,fyPos,cpPos,fCapPos,capSlope,fCapRefPos,LP,dLimPos,fLimPos,snHor,resSnHor,elstk,Resfac)
		327	snCalc();
		328
		329	if((abs(dmax-dyieldPos)>=1.0e-10)&&(abs(sn-dyieldPos)<=1.0e-10)) {
		330	sn=dyieldPos-1.0e-9;
		331	}
		332	}
		333	}
		334
		335	if (sn>dmax) {
		336	dmax = sn;
		337	}
		338
		339	if ((iNoFneg==1)&&(dP<=dNewLoadNeg)) {
		340	sn = dNewLoadNeg - 1.0e-6;
		341	resSn = 0;
		342	}
		343
		344	}
		345	// LOADING ----------------------------------------------------------
		346	// Push envelope
		347	// Compute force and stiffness for this increment
		348	if ((iNoFneg==1)&&(iNoFpos==1)&&(d<fracDispPos)) {
		349	f = Resfac*fyieldPos;
		350	ek = 1.0e-7;
		351	jcode = 8;
		352	} else if((iNoFneg==1)&&(iNoFpos==1)&&(d>=fracDispPos)\|\|flagControlResponse==1) {
		353	f = 1.0e-10;
		354	ek = 1.0e-7;
		355	jcode = 9;
		356	flagstopdeg = 1;
		357	} else if ((iNoFpos==1)&&(d>sn)&&(d<fracDispPos)) {
		358	f = Resfac*fyieldPos;
		359	ek = 1.0e-7;
		360	jcode = 8;
		361	} else if ((iNoFpos==1)&&(d>sn)&&(d>=fracDispPos)){
		362	f = 1.0e-10;
		363	ek =1.0e-7;
		364	jcode = 9;
		365	flagstopdeg = 1;
		366	} else if (d>=fracDispPos \|\| flagControlResponse == 1){
		367	f = 1.0e-10;
		368	ek =1.0e-7;
		369	jcode = 9;
		370	flagstopdeg = 1;
		371	flagControlResponse = 1; // Switches the response to zero strength
		372	} else if ((iNoFpos==1)&&(d<sn)) {
		373	ek = ekunload;
		374	f = fP+ek*deltaD;
		375	jcode = 2;
		376	} else if ((iNoFneg==1)&&(d<dNewLoadNeg)){
		377	f = 0;
		378	ek = 1.0e-7;
		379	jcode = 8;
		380	} else if (d>dmax) {
		381	// call envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac)
		382	f=0.0;
		383	envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac);
		384	dmax = d;
		385	fmax = f;
		386	jcode = 5;
		387	if(ek==elstk*capSlope) {
		388	jcode = 7;
		389	}
		390	if(ek<=1.0e-7) {
		391	jcode = 8;
		392	iCapPos = 1;
		393	}
		394
		395	// c COMPUTE MAXIMUM POSSIBLE DISPLACEMENT
		396	// call boundPos (dBoundPos,fCapRefPos,capSlope,fyNeg,dyNeg,alphaNeg,fCapPos,cpPos,elstk)
		397	dBoundPos=boundPos();
		398	if((d>dBoundPos)\|\|(ek==1.0e-7)) {
		399	iNoFpos = 1;
		400	}
		401
		402	} else if (abs(sn)>1.0e-10) {
		403
		404
		405
		406	if (LP==0) {
		407
		408
		409	if (cpPos<=dyPos) {
		410	// call interPoint(xDevPos1,yDevPos1,sn,resSn,ekhardPos,cpPos,fCapPos,capSlope*elstk)
		411	interPoint(xDevPos1,yDevPos1,sn,resSn,ekhardPos,cpPos,fCapPos,capSlope*elstk);
		412	// call interPoint(xDevPos2,yDevPos2,sn,resSn,ekunload,cpPos,fCapPos,capSlope*elstk)
		413	interPoint(xDevPos2,yDevPos2,sn,resSn,ekunload,cpPos,fCapPos,capSlope*elstk);
		414	if (xDevPos1>xDevPos2) {
		415	xDevPos = xDevPos1;
		416	yDevPos = yDevPos1;
		417	} else {
		418	xDevPos = xDevPos2;
		419	yDevPos = yDevPos2;
		420	}
		421
		422	if ((d<=sn)&&(d<=xDevPos)) {
		423	ek = ekunload;
		424	f = fP+ek*deltaD;
		425	jcode = 2;
		426	} else if ((d>sn)&&(d<=xDevPos)) {
		427	ek = elstk*alphaPos;
		428	f2 = resSn+ek*(d-sn);
		429	f1 = fP+ekunload*deltaD ;
		430	//f = min(f1,f2);
		431	if (f1<f2)
		432	{
		433	f=f1;
		434	}
		435	else
		436	{
		437	f=f2;
		438	}
		439
		440	jcode = 5;
		441
		442	} else if (d>xDevPos) {
		443	ek = capSlope*elstk;
		444	f2 = yDevPos+ek*(d-xDevPos);
		445	f1 = fP+ekunload*deltaD ;
		446	//f = min(f1,f2);
		447	if (f1<f2)
		448	{
		449	f=f1;
		450	}
		451	else
		452	{
		453	f=f2;
		454	}
		455	if (ek!=ekunload) {
		456	// call envHitsZero(f,fP,ek)
		457	envHitsZero(f);
		458	}
		459	iCapPos = 1;
		460	jcode = 7;
		461	} else {
		462	}
		463
		464	} else if (cpPos > dyPos) {
		465	interPoint(xDevPos1,yDevPos1,sn,resSn,ekhardPos,cpPos,fCapPos,capSlope*elstk);
		466	if(d<=sn && d<=xDevPos1) { // Unloading in positive loading direction
		467	ek = ekunload;
		468	f = fP+ek*deltaD;
		469	jcode = 2;
		470	} else if ((d>sn)&&(d<=cpPos)&& (d<=xDevPos1)) {
		471	ek = elstk*alphaPos;
		472	f2 = resSn+ek*(d-sn);
		473	f1 = fP+ekunload*deltaD;
		474	//f = min(f1,f2);
		475	if (f1<f2)
		476	{
		477	f=f1;
		478	}
		479	else
		480	{
		481	f=f2;
		482	}
		483	if (abs(f-f1)<1.0e-10) {
		484	ek=ekunload;
		485	}
		486	jcode = 5;
		487	} else if((d>cpPos)&&(d<fracDispPos)) {
		488	ek = capSlope*elstk;
		489	f2 = fCapPos+ek*(d-cpPos);
		490	f1 = fP+ekunload*deltaD;
		491	//f = min(f1,f2);
		492	if (f1<f2)
		493	{
		494	f=f1;
		495	}
		496	else
		497	{
		498	f=f2;
		499	}
		500	if (abs(f-f1)<1.0e-10) {
		501	ek=ekunload;
		502	}
		503	if (ek!=ekunload) {
		504	// call envHitsZero(f,fP,ek)
		505	envHitsZero(f);
		506	}
		507	iCapPos = 1;
		508	jcode = 7;
		509	// c added by Dimitrios to simulate ductile fracture
		510	} else if(d>=fracDispPos){
		511	f=1.0e-10;
		512	ek = -1.0e-7;
		513	jcode = 9;
		514	flagstopdeg = 1;
		515	}
		516	// Added by Dimitrios to stay on the residual path when dresidual <d< dfracture
		517	if(d < fracDispPos && d > cpPos+(ResfacfyieldPos-fCapPos)/(capSlopeelstk)) {
		518	f = Resfac*fyieldPos;
		519	ek = -1.0e-7;
		520	jcode = 9;
		521	}
		522	if( flagControlResponse == 1) { // Response has to be zero since post capping regions hit zero
		523	f = 1.0e-10;
		524	jcode = 10; // this is for DRAIN output path tracking
		525	}
		526	}
		527
		528	// c IF LP IS EQUAL TO 1
		529	} else if (LP==1) {
		530	if(d<=sn) {
		531	ek = ekunload;
		532	f = fP+ek*deltaD;
		533	jcode = 2;
		534	} else if (((d>sn)&&(sn==snEnv)&&(d<=snHor))\|\|((iNoFneg==1)&&(d>sn)&&(d<snHor))) {
		535	ek = elstk*alphaPos;
		536	f2 = resSn+ek*(d-sn);
		537	f1 = fP+ekunload*deltaD;
		538	//f = min(f1,f2);
		539	if (f1<f2)
		540	{
		541	f=f1;
		542	}
		543	else
		544	{
		545	f=f2;
		546	}
		547	if (abs(f-f1)<1.0e-10) {
		548	ek=ekunload;
		549	}
		550	jcode = 5;
		551
		552	} else {
		553	ek = 0;
		554	f1 = fP+ekunload*deltaD;
		555	f2 = fLimPos;
		556	//f = min(f1,f2);
		557	if (f1<f2)
		558	{
		559	f=f1;
		560	}
		561	else
		562	{
		563	f=f2;
		564	}
		565	if (abs(f-f1)<1.0e-10) {
		566	ek=ekunload;
		567	}
		568	jcode = 8;
		569	}
		570	} else {
		571	}
		572	// c Elastic
		573	} else {
		574	if (d>0.0) {
		575	// call envelPosCap2 (fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac)
		576	// [d,f,ek]=envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,0.0,ek,elstk,fyieldPos,Resfac);
		577	f=0.0;
		578	envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac);
		579	} else {
		580	// call envelNegCap2 (fyNeg,alphaNeg,capSlope,cpNeg,d,f,ek,elstk,fyieldNeg,Resfac)
		581	f=0.0;
		582	envelNegCap2(fyNeg,alphaNeg,capSlopeNeg,cpNeg,d,f,ek,elstk,fyieldNeg,ResfacNeg);
		583	}
		584	jcode = 0;
		585	if(ek==elstk*alphaPos) {
		586	jcode = 1;
		587	}
		588	if(ek==elstk*capSlope) {
		589	jcode = 7;
		590	}
		591	if(ek<=1.0e-7) {
		592	jcode = 8;
		593	}
		594	}
		595
		596	//// c IF D E L T A < 0 - - - - - - - - - - - - - - - - - - - - - - - -
		597
		598	} else {
		599
		600	if (iNoFneg==1) {
		601	// call interPoint(dNewLoadNeg,fNewLoadNeg,dyPos,fyPos,ekhardPos,0.d0,0.d0,0.d0)
		602	interPoint(dNewLoadNeg,fNewLoadNeg,dyPos,fyPos,ekhardPos,0.0,0.0,0.0);
		603	}
		604
		605
		606	// c If there is a corner changing delta from positive to negative ---
		607
		608	if (kon==1) {
		609	kon = 2;
		610
		611	dlstPos = dP;
		612	flstPos = fP;
		613
		614	if (dlstPos>=dmax) {
		615	fLimPos = flstPos;
		616	dLimPos = dlstPos;
		617	}
		618
		619	if (resSp<0) {
		620	RSE = 0.5fPfP/ekunload;
		621	} else {
		622	RSE=0.5(fP+resSn)(dP-sp)+0.5resSpresSp/(elstk*alphaNeg);
		623	}
		624
		625	a2 = Enrgtk-(Enrgtot-RSE);
		626
		627
		628	if((a2<=0.0)&&(Enrgtk!=0.0)) {
		629	flgstop=1;
		630	}
		631	//// Update the Unloading Stiffness deterioration
		632	if((ecapk!=0.0)&&(d>sn)&&(flagstopdeg!=1)&&(flagdeg!=1)) {
		633	betak = pow(((Enrgc-RSE)/(Enrgtk-(Enrgtot-RSE))),ck);
		634	if(((Enrgtot-RSE)>=Enrgtk)\|\|(betak>=1.0)) {
		635	betak = 1.0;
		636	}
		637	// If Post caping slopes have not been flat due to residual update stiffness deterioration
		638	if(flagdeg !=1 \|\| flagstopdeg!=1) {
		639	ekunload = ekexcurs*(1-betak);
		640	} else { // Keep the same unloading stiffness
		641	ekunload = ekexcurs;
		642	}
		643	if(ekunload<=0.1*elstk) {
		644	ekunload = 0.1*elstk;
		645	}
		646	}
		647
		648
		649	//// sp CALCULATION----------------------------------------------------
		650
		651	if((dmin<dyieldNeg)\|\|(dmax>dyieldPos)) {
		652	if((dP>sn)\|\|(((dP<sn)&&(ekP==ekunload)))) {
		653
		654	// call spCalc(sp,resSp,spEnv,resSpEnv,dP,fP,ekunload,alphaNeg,dyNeg,fyNeg,cpNeg,fCapNeg,capSlope,fCapRefNeg,LN,dLimNeg,fLimNeg,spHor,resSpHor,elstk,Resfac)
		655	spCalc();
		656	if((abs(dmin-dyieldNeg)>=1.0e-10)&&(abs(sp-dyieldNeg)<=1.0e-10)) {
		657	sp=dyieldNeg-1.0e-9;
		658	}
		659	}
		660	}
		661
		662	if (sp<dmin) {
		663	dmin = sp;
		664	}
		665
		666	if ((iNoFpos==1)&&(dP>=dNewLoadPos)) {
		667	sp = dNewLoadPos + 1.0e-6;
		668	resSp = 0;
		669	}
		670	}
		671
		672	//// UNLOADING
		673	// c Push envelope
		674
		675	if ((iNoFneg==1)&&(iNoFpos==1)&&(d>fracDispNeg)) {
		676	f = ResfacNeg*fyieldNeg;
		677	ek = 1.0e-7;
		678	jcode = 8;
		679	} else if ((iNoFneg==1)&&(iNoFpos==1)&&(d<=fracDispNeg)\|\|flagControlResponse==1){
		680	f = -1.0e-10;
		681	ek =1.0e-7;
		682	jcode = 9;
		683	flagstopdeg = 1;
		684	} else if ((iNoFneg==1)&&(d<sp)&&(d>fracDispNeg)) {
		685	f = ResfacNeg*fyieldNeg;
		686	ek = 1.0e-7;
		687	jcode = 8;
		688	} else if ((iNoFneg==1)&&(d<sp)&&(d<=fracDispNeg)) {
		689	f = -1.0e-10;
		690	ek = 1.0e-7;
		691	jcode = 9;
		692	flagstopdeg = 1;
		693	} else if (d<=fracDispNeg \|\| flagControlResponse ==1){
		694	f = -1.0e-10;
		695	ek = 1.0e-7;
		696	jcode = 9;
		697	flagstopdeg = 1;
		698	flagControlResponse = 1; // To control response after I exceed fracture
		699	} else if ((iNoFneg==1)&&(d>=sp)) {
		700	ek = ekunload;
		701	f = fP+ek*deltaD;
		702	jcode = 2;
		703	} else if ((iNoFpos==1)&&(d>dNewLoadPos)) {
		704	f = 0;
		705	ek = 1.0e-7;
		706	jcode = 8;
		707	} else if (d<dmin) {
		708	// call envelNegCap2(fyNeg,alphaNeg,capSlope,cpNeg,d,f,ek,elstk,fyieldNeg,Resfac)
		709	f=0.0;
		710	envelNegCap2(fyNeg,alphaNeg,capSlopeNeg,cpNeg,d,f,ek,elstk,fyieldNeg,ResfacNeg);
		711	dmin = d;
		712	fmin = f;
		713
		714	jcode = 5;
		715	if(ek==elstk*capSlopeNeg) {
		716	jcode = 7;
		717	}
		718	if(ek<=1.0e-7) {
		719	jcode = 8;
		720	iCapNeg = 1;
		721	}
		722
		723	//// c COMPUTE MINIMUM POSSIBLE DISPLACEMENT
		724	// call boundNeg (dBoundNeg,fCapRefNeg,capSlope,fyPos,dyPos,alphaPos,fCapNeg,cpNeg,elstk)
		725	dBoundNeg=boundNeg();
		726	if((d<dBoundNeg)\|\|(ek==1.0e-7)) {
		727	iNoFneg = 1;
		728	}
		729
		730	} else if (abs(sp)>1.0e-10){
		731
		732	// c If LN is equal to zero
		733	if (LN==0) {
		734
		735	if (cpNeg>=dyNeg) {
		736	// call interPoint(xDevNeg1,yDevNeg1,sp,resSp,elstkalphaNeg,cpNeg,fCapNeg,capSlopeelstk)
		737	interPoint(xDevNeg1,yDevNeg1,sp,resSp,elstkalphaNeg,cpNeg,fCapNeg,capSlopeNegelstk);
		738	// call interPoint(xDevNeg2,yDevNeg2,sp,resSp,ekunload,cpNeg,fCapNeg,capSlope*elstk)
		739	interPoint(xDevNeg2,yDevNeg2,sp,resSp,ekunload,cpNeg,fCapNeg,capSlopeNeg*elstk);
		740	if (xDevNeg1<xDevNeg2) {
		741	xDevNeg = xDevNeg1;
		742	yDevNeg = yDevNeg1;
		743	} else {
		744	xDevNeg = xDevNeg2;
		745	yDevNeg = yDevNeg2;
		746	}
		747	//
		748	if ((d>=sp)&&(d>=xDevNeg)) {
		749	ek = ekunload;
		750	f = fP+ek*deltaD;
		751	jcode = 2;
		752	} else if ((d<sp)&&(d>=xDevNeg)) {
		753	ek = elstk*alphaNeg;
		754	f2 = resSp+ek*(d-sp);
		755	f1 = fP+ekunload*deltaD;
		756	//f = max(f1,f2);
		757	if (f1>f2)
		758	{
		759	f=f1;
		760	}
		761	else
		762	{
		763	f=f2;
		764	}
		765	if (abs(f-f1)<1.0e-10) {
		766	ek=ekunload;
		767	}
		768	jcode = 5;
		769	} else if (d<xDevNeg) {
		770	ek = capSlopeNeg*elstk;
		771	f2 = yDevNeg+ek*(d-xDevNeg);
		772	f1 = fP+ekunload*deltaD;
		773	//f = max(f1,f2);
		774	if (f1>f2)
		775	{
		776	f=f1;
		777	}
		778	else
		779	{
		780	f=f2;
		781	}
		782	if (abs(f-f1)<1.0e-10) {
		783	ek=ekunload;
		784	}
		785	if (ek!=ekunload) {
		786	// call envHitsZero(f,fP,ek)
		787	envHitsZero(f);
		788	}
		789	iCapNeg = 1;
		790	jcode = 7;
		791	// } else
		792	}
		793
		794
		795	} else if (cpNeg<dyNeg) {
		796	interPoint(xDevNeg1,yDevNeg1,sp,resSp,elstkalphaNeg,cpNeg,fCapNeg,capSlopeNegelstk);
		797	if ( d>=sp && d>=xDevNeg1) {
		798	ek = ekunload;
		799	f = fP+ek*deltaD;
		800	jcode = 2;
		801	} else if((d<sp)&&(d>=cpNeg)&& (d>=xDevNeg1)) {
		802	ek = elstk*alphaNeg;
		803	f2 = resSp+ek*(d-sp);
		804	f1 = fP+ekunload*deltaD;
		805	//f=max(f1,f2);
		806	if (f1>f2)
		807	{
		808	f=f1;
		809	}
		810	else
		811	{
		812	f=f2;
		813	}
		814	if (abs(f-f1)<1.0e-10) {
		815	ek=ekunload;
		816	}
		817	jcode = 5;
		818	} else if((d<cpNeg)&&(d>fracDispNeg)) {
		819	ek = capSlopeNeg*elstk;
		820	f2 = fCapNeg+ek*(d-cpNeg);
		821	f1 = fP+ekunload*deltaD;
		822	//f = max(f1,f2);
		823	if (f1>f2)
		824	{
		825	f=f1;
		826	}
		827	else
		828	{
		829	f=f2;
		830	}
		831	if (abs(f-f1)<1.0e-10) {
		832	ek=ekunload;
		833	}
		834	if (ek!=ekunload) {
		835	// call envHitsZero(f,fP,ek)
		836	envHitsZero(f);
		837	}
		838	iCapNeg = 1;
		839	jcode = 7;
		840	// Added by Dimitris to model ductile tearing. Once theta_u(fracDispNeg) is exceeded Strength drops to zero and stays
		841	} else if(d<=fracDispNeg \|\| flagControlResponse == 1){
		842	ek = -1.0e-7;
		843	f = 1.0e-10;
		844	jcode = 9;
		845	flagstopdeg = 1;
		846	flagControlResponse = 1; // To dictate the response after passing fracture
		847	}
		848	// Added by Dimitrios to stay on the residual path when dresidual <d- < dfracture
		849	if(d > fracDispNeg && d < cpNeg+(ResfacNegfyieldNeg-fCapNeg)/(capSlopeNegelstk)) {
		850	f = ResfacNeg*fyieldNeg;
		851	ek = -1.0e-7;
		852	jcode = 9;
		853	}
		854	if( flagControlResponse == 1) { // Response has to be zero since post capping regions hit zero
		855	f = 1.0e-10;
		856	jcode = 10; // this is for DRAIN output path tracking
		857	}
		858	}
		859
		860	} else if (LN==1){
		861	if(d>=sp){
		862	ek = ekunload;
		863	f = fP+ek*deltaD;
		864	jcode = 2;
		865	} else if (((d<sp)&&(sp==spEnv)&&(d>spHor))\|\|((iNoFpos==1)&&(d<sp)&&(d>spHor))) {
		866	ek = elstk*alphaNeg;
		867	f2 = resSp+ek*(d-sp);
		868	f1 = fP+ekunload*deltaD;
		869	//f = max(f1,f2);
		870	if (f1>f2)
		871	{
		872	f=f1;
		873	}
		874	else
		875	{
		876	f=f2;
		877	}
		878	if (abs(f-f1)<1.0e-10) {
		879	ek=ekunload;
		880	}
		881	jcode = 5;
		882	} else {
		883	ek = 1.0e-7;
		884	f1 = fP+ekunload*deltaD ;
		885	f2 = fLimNeg;
		886	//f = max(f1,f2);
		887	if (f1>f2)
		888	{
		889	f=f1;
		890	}
		891	else
		892	{
		893	f=f2;
		894	}
		895	if (abs(f-f1)<1.0e-10) {
		896	ek=ekunload;
		897	}
		898	jcode = 8;
		899	}
		900	} else {
		901	}
		902	} else {
		903	if (d>0.0) {
		904	// call envelPosCap2 (fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac)
		905	//[d,f,ek]=envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,0.0,ek,elstk,fyieldPos,Resfac);
		906	f=0.0;
		907	envelPosCap2(fyPos,alphaPos,capSlope,cpPos,d,f,ek,elstk,fyieldPos,Resfac);
		908	} else {
		909	// call envelNegCap2 (fyNeg,alphaNeg,capSlope,cpNeg,d,f,ek,elstk,fyieldNeg,Resfac)
		910	envelNegCap2(fyNeg,alphaNeg,capSlopeNeg,cpNeg,d,f,ek,elstk,fyieldNeg,ResfacNeg);
		911	}
		912	jcode = 0;
		913	if(ek==elstk*alphaNeg) {
		914	jcode = 1;
		915	}
		916	if(ek==elstk*capSlopeNeg) {
		917	jcode = 7;
		918	}
		919	if(ek<=1.0e-7) {
		920	jcode = 8;
		921	}
		922	}
		923
		924	}
		925	// c }S BIG LOOP ****************************************************
		926	////
		927	//
		928	dn = d;
		929	fn = f;
		930
		931	Enrgi = 0.5(f+fP)deltaD;
		932
		933	Enrgc = Enrgc + Enrgi;
		934	Enrgtot = Enrgtot + Enrgi;
		935
		936	RSE = 0.5ff/ekunload;
		937
		938	//// Flag to deteriorate parameters on the opposite side of the loop --
		939
		940	if((f*fP<0.0)&&(interup==0)) {
		941	if(((fP>0.0)&&(dmax>dyieldPos))\|\|((fP<0.0)&&(dmin<dyieldNeg))) {
		942	flagdeg = 1;
		943	interup=1;
		944	}
		945	}
		946
		947
		948	//// energy CALCULATIONS ---------------------------------------------
		949
		950	if((flagstopdeg==0)&&(flagdeg==1)) {
		951	iDeg = iDeg +1;
		952
		953	if((Enrgtot>=Enrgts)&&(Enrgts!=0.0)) {
		954	betas = 1.0;
		955	} else if((Enrgtot>=Enrgtd)&&(Enrgtd!=0.0)) {
		956	betad = 1.0;
		957	} else {
		958	if(ecaps!=0.0) {
		959	betas = pow((Enrgc/(Enrgts-Enrgtot)),cs);
		960	}
		961	if(ecapd!=0.0) {
		962	betad = pow((Enrgc/(Enrgtd-Enrgtot)),cd);
		963	}
		964
		965	if(abs(betas)>=1.0) {
		966	betas = 1.0;
		967	}
		968	if(abs(betad)>=1.0) {
		969	betad = 1.0;
		970	}
		971	}
		972	//// Initialize energy of the cycle and Kstif for next loop -------
		973
		974	Enrgc = 0.0;
		975	ekexcurs = ekunload;
		976
		977	//// Check the horizontal limit for subsequent cycles -------------
		978	if ((dmin<cpNeg)&&(LN==1)) {
		979	iCapNeg = 1;
		980	}
		981	if ((dmax>cpPos)&&(LP==1)) {
		982	iCapPos = 1;
		983	}
		984
		985
		986	//// Deteriorate parameters for the next half cycle
		987	if(deltaD<0.0){
		988
		989	if(flagstopdeg == 0){
		990	fyNeg = fyNeg(1-betasDNeg);
		991	alphaNeg=alphaNeg(1-betasDNeg);
		992	fCapRefNeg=fCapRefNeg(1-betadDNeg);
		993	}else{
		994	fyNeg = fyNeg;
		995	alphaNeg=alphaNeg;
		996	fCapRefNeg=fCapRefNeg;
		997	}
		998	// When we reach post capping slope goes to zero due to residual
		999	if(fyNeg>=ResfacNeg*fyieldNeg) { // If strength drops below residual
		1000	fyNeg = ResfacNeg*fyieldNeg;
		1001	alphaNeg = 10^(-4);
		1002	fCapRefNeg = fyNeg;
		1003	capSlopeNeg = -pow(10.0,-6);
		1004	flagstopdeg = 1;
		1005	} else { //% Keep updating the post capping slope
		1006	capSlopeNeg = capSlopeOrigNeg(1-abs((ResfacNegfyieldNeg)/fyNeg));
		1007	if(capSlopeNeg >=0){
		1008	capSlopeNeg = -pow(10.0,-6);
		1009	}
		1010	}
		1011
		1012	dyNeg = fyNeg/elstk;
		1013	ekhardNeg=alphaNeg*elstk;
		1014
		1015	dCap1Neg=fCapRefNeg/(elstk-capSlopeOrigNeg*elstk);
		1016	dCap2Neg=(fCapRefNeg+ekhardNegdyNeg-fyNeg)/(ekhardNeg-capSlopeOrigNegelstk);
		1017	//cpNeg=min(dCap1Neg,dCap2Neg);
		1018	if (dCap1Neg<dCap2Neg)
		1019	{
		1020	cpNeg=dCap1Neg;
		1021	}
		1022	else
		1023	{
		1024	cpNeg=dCap2Neg;
		1025	}
		1026
		1027	fCapNeg = fCapRefNeg + capSlopeOrigNegelstkcpNeg;
		1028
		1029	// call envelNegCap2 (fyNeg,alphaNeg,capSlope,cpNeg,dLimNeg,fLimNeg,ekt,elstk,fyieldNeg,Resfac)
		1030	envelNegCap2(fyNeg,alphaNeg,capSlopeNeg,cpNeg,dLimNeg,fLimNeg,ekt,elstk,fyieldNeg,ResfacNeg);
		1031	// call spCalc(sp,resSp,spEnv,resSpEnv,dP,fP,ekunload,alphaNeg,dyNeg,fyNeg,cpNeg,fCapNeg,capSlope,fCapRefNeg,LN,dLimNeg,fLimNeg,spHor,resSpHor,elstk,Resfac)
		1032	spCalc();
		1033	// c In case the degradation point moves from the negative to pos. side
		1034	if(resSp>f) {
		1035	d = sp;
		1036	f = resSp;
		1037	ek = ekhardNeg;
		1038	}
		1039	} else {
		1040	if(flagstopdeg == 0){
		1041	fyPos = fyPos(1-betasDPlus);
		1042	alphaPos=alphaPos(1-betasDPlus);
		1043	fCapRefPos=fCapRefPos(1-betadDPlus);
		1044	} else {
		1045	fyPos = fyPos;
		1046	alphaPos=alphaPos;
		1047	fCapRefPos=fCapRefPos;
		1048	}
		1049
		1050	// %If post capping slope goes to zero due to residual:
		1051	if(fyPos <= Resfac*fyieldPos) { //% If yield Strength Pos drops below residual
		1052	fyPos = Resfac*fyieldPos;
		1053	alphaPos = pow(10.0,-4);
		1054	fCapRefPos = fyPos;
		1055	capSlope = -pow(10.0,-6);
		1056	flagstopdeg = 1;
		1057	} else { //% keep updating
		1058	capSlope = capSlopeOrig(1-abs((ResfacfyieldPos)/fyPos));
		1059	if(capSlope >=0) {
		1060	capSlope = -pow(10.0,-6);
		1061	}
		1062	}
		1063	dyPos = fyPos/elstk;
		1064	ekhardPos=alphaPos*elstk;
		1065
		1066	dCap1Pos=fCapRefPos/(elstk-capSlopeOrig*elstk);
		1067	dCap2Pos=(fCapRefPos+ekhardPosdyPos-fyPos)/(ekhardPos-capSlopeOrigelstk);
		1068	//cpPos=max(dCap1Pos,dCap2Pos);
		1069	if(dCap1Pos>dCap2Pos)
		1070	{
		1071	cpPos=dCap1Pos;
		1072	}
		1073	else
		1074	{
		1075	cpPos=dCap2Pos;
		1076	}
		1077	fCapPos = fCapRefPos + capSlopeOrigelstkcpPos;
		1078
		1079	//call envelPosCap2 (fyPos,alphaPos,capSlope,cpPos,dLimPos,fLimPos,ekt,elstk,fyieldPos,Resfac)
		1080	envelPosCap2(fyPos,alphaPos,capSlope,cpPos,dLimPos,fLimPos,ekt,elstk,fyieldPos,Resfac);
		1081
		1082	//call snCalc(sn,resSn,snEnv,resSnEnv,dP,fP,ekunload,alphaPos,dyPos,fyPos,cpPos,fCapPos,capSlope,fCapRefPos,LP,dLimPos,fLimPos,snHor,resSnHor,elstk,Resfac)
		1083	snCalc();
		1084	// c In case the degradation point moves from the pos. to neg. side
		1085	if(resSn<f) {
		1086	d = sn;
		1087	f = resSn;
		1088	ek = ekhardPos;
		1089
		1090	}
		1091	}
		1092	}
		1093
		1094	// c Check the horizontal limit in case that dBound is reached after first neg slope
		1095	if ((d<0)&&(abs(ek)<=1.0e-7)) {
		1096	LN = 1;
		1097	}
		1098	if ((d>0)&&(abs(ek)<=1.0e-7)) {
		1099	LP = 1;
		1100	}
		1101
		1102
		1103	// c Calculation of static resisting force vector for the element -----
		1104
		1105	double RestoringForce=f;
		1106
		1107
		1108
		1109
		1110
		1111	// c Update envelope values --------------------------------------------
		1112
		1113	f1 = ftot;
		1114	v1 = vtot;
		1115	d1 = dtot;
		1116
		1117	ftot = fn;
		1118	dtot = dn;
		1119	vtot = dn;
		1120
		1121
		1122
		1123
		1124
		1125	//// c Change in Energies -----------------------------------------------
		1126
		1127	// c Hysteretic energy
		1128	ener = Enrgi;
		1129
		1130	// c Updating parameters for next cycle ---------------------------------
		1131	stif = ek;
		1132	ekP = ek;
		1133	fP = f;
		1134	dP = d;
		1135
		1136	if(jcode!=kcode) {
		1137	kst = 1;
		1138	}
		1139	kcode = jcode;
		1140
		1141	//priority of logical operators
		1142	if (interup==1&&(ek==elstkalphaPos) \|\|(ek==elstkalphaNeg)\|\|(ek==capSlopeelstk) \|\|(ek==capSlopeNegelstk)) {
		1143	interup = 0;
		1144	}
		1145
		1146
		1147
		1148	//*****************************} of matlab code*****************************************************
		1149
		1150	//Define force and tangent
		1151	Force=RestoringForce;
		1152	Tangent=stif;
		1153	return 0;
		1154	}
		1155	double
		1156	Bilin::getStress(void)
		1157	{
		1158	return (Force);
		1159	}
		1160
		1161	double
		1162	Bilin::getTangent(void)
		1163	{
		1164
		1165	return (Tangent);
		1166	}
		1167
		1168	double
		1169	Bilin::getInitialTangent(void)
		1170	{
		1171	return (Ke);
		1172	}
		1173
		1174
		1175
		1176
		1177	double
		1178	Bilin::getStrain(void)
		1179	{
		1180	return (U);
		1181	}
		1182
		1183
		1184
		1185	int
		1186	Bilin::commitState(void)
		1187	{
		1188	//commit trial variables
		1189	//CU=TU; //displacement
		1190	CU=U;
		1191	CForce=Force;
		1192	CTangent=Tangent;
		1193
		1194	CdNewLoadPos=dNewLoadPos;
		1195	CdNewLoadNeg=dNewLoadNeg;
		1196	Cflgstop=flgstop;
		1197	Cflagdeg=flagdeg ;
		1198	Cflagstopdeg=flagstopdeg;
		1199	Cekt=ekt;
		1200	Cinterup=interup;
		1201	Ckcode=kcode;
		1202	Ckon=kon;
		1203	CiCapNeg=iCapNeg;
		1204	CiNoFneg=iNoFneg;
		1205	CiNoFpos=iNoFpos;
		1206	CiCapPos=iCapPos;
		1207	CiDeg=iDeg;
		1208	CLP=CLP;
		1209	CLN=LN;
		1210	CcapSlope=capSlope;
		1211	CcapDispPos=capDispPos;
		1212	CcapDispNeg=capDispNeg;
		1213	Celstk=elstk;
		1214	CfyieldPos=fyieldPos;
		1215	CfyieldNeg=fyieldNeg;
		1216	Calpha=alpha;
		1217	Cecaps=ecaps;
		1218	Cecapk=ecapk;
		1219	Cecapd=ecapd;
		1220	Ccs=cs;
		1221	Cck=ck;
		1222	Ccd=cd;
		1223	Cdmax=dmax;
		1224	Cdmin=dmin;
		1225	CEnrgtot=Enrgtot;
		1226	CEnrgc=Enrgc;
		1227	CfyPos=fyPos;
		1228	CfLimNeg=fLimNeg;
		1229	CfyNeg=fyNeg;
		1230	CekP=ekP;
		1231	Cekunload=ekunload;
		1232	Csp=sp;
		1233	Csn=sn;
		1234	CdP=dP;
		1235	CfP=fP;
		1236	Cek=ek;
		1237	Cstif=stif;
		1238	CdLimPos=dLimPos;
		1239	CdLimNeg=dLimNeg;
		1240	Cvtot=vtot;
		1241	Cftot=ftot;
		1242	Cdtot=dtot;
		1243	Cdn=dn;
		1244	CcpPos=cpPos;
		1245	CcpNeg=cpNeg;
		1246	CfLimPos=fLimPos;
		1247	CdlstPos=dlstPos;
		1248	CflstPos=flstPos;
		1249	CdlstNeg=dlstNeg;
		1250	CflstNeg=flstNeg;
		1251	Cekexcurs=ekexcurs;
		1252	CRSE=RSE;
		1253	CfPeakPos=fPeakPos;
		1254	CfPeakNeg=fPeakNeg;
		1255	CdCap1Pos=dCap1Pos;
		1256	CdCap2Pos=dCap2Pos;
		1257	CdCap1Neg=dCap1Neg;
		1258	CdCap2Neg=dCap2Neg;
		1259	CalphaNeg=alphaNeg;
		1260	CalphaPos=alphaPos;
		1261	CekhardNeg=ekhardNeg;
		1262	CekhardPos=ekhardPos;
		1263	CfCapRefPos=fCapRefPos;
		1264	CfCapRefNeg=fCapRefNeg;
		1265	CEnrgts=Enrgts;
		1266	CEnrgtk=Enrgtk;
		1267	CEnrgtd=Enrgtd;
		1268	CdyPos=dyPos;
		1269	CdyNeg=dyNeg;
		1270	CdyieldPos=dyieldPos;
		1271	CdyieldNeg=dyieldNeg;
		1272	CresSnHor=resSnHor;
		1273	Cfmax=fmax;
		1274	Cfmin=fmin;
		1275	CresSp=resSp;
		1276	CresSn=resSn;
		1277	CfCapPos=fCapPos;
		1278	CfCapNeg=fCapNeg;
		1279	CsnHor=snHor;
		1280	CspHor=spHor;
		1281	CresSpHor=resSpHor;
		1282	CsnEnv=snEnv;
		1283	CresSnEnv=resSnEnv;
		1284	CspEnv=spEnv;
		1285	CresSpEnv=resSpEnv;
		1286	CResfac=Resfac;
		1287	CcapSlopeOrig=capSlopeOrig;
		1288	CfracDispPos=fracDispPos;
		1289	CfracDispNeg=fracDispNeg;
		1290	CDPlus=DPlus;
		1291	CDNeg=DNeg;
		1292	CalphaN=alphaN;
		1293	Cecapa=ecapa;
		1294	Cca=ca;
		1295	CResfacNeg=ResfacNeg;
		1296	CmyFcapping=myFcapping;
		1297	CmyFcappingNeg=myFcappingNeg;
		1298	CcapSlopeNeg=capSlopeNeg;
		1299	CflagControlResponse=flagControlResponse;
		1300	CcapSlopeOrigNeg=capSlopeOrigNeg;
		1301	CUprev=Uprev;
		1302	return 0;
		1303	}
		1304
		1305	int
		1306	Bilin::revertToLastCommit(void)
		1307	{
		1308	//the opposite of commit trial history variables
		1309	U=CU;
		1310	Force=CForce;
		1311	Tangent=CTangent;
		1312
		1313	dNewLoadPos=CdNewLoadPos;
		1314	dNewLoadNeg=CdNewLoadNeg;
		1315	flgstop=Cflgstop;
		1316	flagdeg=Cflagdeg ;
		1317	flagstopdeg=Cflagstopdeg;
		1318	ekt=Cekt;
		1319	interup=Cinterup;
		1320	kcode=Ckcode;
		1321	kon=Ckon;
		1322	iCapNeg=CiCapNeg;
		1323	iNoFneg=CiNoFneg;
		1324	iNoFpos=CiNoFpos;
		1325	iCapPos=CiCapPos;
		1326	iDeg=CiDeg;
		1327	LP=CLP;
		1328	LN=CLN;
		1329	capSlope=CcapSlope;
		1330	capDispPos=CcapDispPos;
		1331	capDispNeg=CcapDispNeg;
		1332	elstk=Celstk;
		1333	fyieldPos=CfyieldPos;
		1334	fyieldNeg=CfyieldNeg;
		1335	alpha=Calpha;
		1336	ecaps=Cecaps;
		1337	ecapk=Cecapk;
		1338	ecapd=Cecapd;
		1339	cs=Ccs;
		1340	ck=Cck;
		1341	cd=Ccd;
		1342	dmax=Cdmax;
		1343	dmin=Cdmin;
		1344	Enrgtot=CEnrgtot;
		1345	Enrgc=CEnrgc;
		1346	fyPos=CfyPos;
		1347	fLimNeg=CfLimNeg;
		1348	fyNeg=CfyNeg;
		1349	ekP=CekP;
		1350	ekunload=Cekunload;
		1351	sp=Csp;
		1352	sn=Csn;
		1353	dP=CdP;
		1354	fP=CfP;
		1355	ek=Cek;
		1356	stif=Cstif;
		1357	dLimPos=CdLimPos;
		1358	dLimNeg=CdLimNeg;
		1359	vtot=Cvtot;
		1360	ftot=Cftot;
		1361	dtot=Cdtot;
		1362	dn=Cdn;
		1363	cpPos=CcpPos;
		1364	cpNeg=CcpNeg;
		1365	fLimPos=CfLimPos;
		1366	dlstPos=CdlstPos;
		1367	flstPos=CflstPos;
		1368	dlstNeg=CdlstNeg;
		1369	flstNeg=CflstNeg;
		1370	ekexcurs=Cekexcurs;
		1371	RSE=CRSE;
		1372	fPeakPos=CfPeakPos;
		1373	fPeakNeg=CfPeakNeg;
		1374	dCap1Pos=CdCap1Pos;
		1375	dCap2Pos=CdCap2Pos;
		1376	dCap1Neg=CdCap1Neg;
		1377	dCap2Neg=CdCap2Neg;
		1378	alphaNeg=CalphaNeg;
		1379	alphaPos=CalphaPos;
		1380	ekhardNeg=CekhardNeg;
		1381	ekhardPos=CekhardPos;
		1382	fCapRefPos=CfCapRefPos;
		1383	fCapRefNeg=CfCapRefNeg;
		1384	Enrgts=CEnrgts;
		1385	Enrgtk=CEnrgtk;
		1386	Enrgtd=CEnrgtd;
		1387	dyPos=CdyPos;
		1388	dyNeg=CdyNeg;
		1389	dyieldPos=CdyieldPos;
		1390	dyieldNeg=CdyieldNeg;
		1391	resSnHor=CresSnHor;
		1392	fmax=Cfmax;
		1393	fmin=Cfmin;
		1394	resSp=CresSp;
		1395	resSn=CresSn;
		1396	fCapPos=CfCapPos;
		1397	fCapNeg=CfCapNeg;
		1398	snHor=CsnHor;
		1399	spHor=CspHor;
		1400	resSpHor=CresSpHor;
		1401	snEnv=CsnEnv;
		1402	resSnEnv=CresSnEnv;
		1403	spEnv=CspEnv;
		1404	resSpEnv=CresSpEnv;
		1405	Resfac=CResfac;
		1406	capSlopeOrig=CcapSlopeOrig;
		1407	fracDispPos=CfracDispPos;
		1408	fracDispNeg=CfracDispNeg;
		1409	DPlus=CDPlus;
		1410	DNeg=CDNeg;
		1411	alphaN=CalphaN;
		1412	ecapa=Cecapa;
		1413	ca=Cca;
		1414	ResfacNeg=CResfacNeg;
		1415	myFcapping=CmyFcapping;
		1416	myFcappingNeg=CmyFcappingNeg;
		1417	capSlopeNeg=CcapSlopeNeg;
		1418	flagControlResponse=CflagControlResponse;
		1419	capSlopeOrigNeg=CcapSlopeOrigNeg;
		1420	Uprev=CUprev;
		1421	return 0;
		1422	}
		1423
		1424	int
		1425	Bilin::revertToStart(void)
		1426	{
		1427	//initially I zero everything
		1428	U=CU=0.0;
		1429	Force=CForce=0.0;
		1430	Tangent=CTangent=0.0;
		1431
		1432	dNewLoadPos=CdNewLoadPos=0.0;
		1433	dNewLoadNeg=CdNewLoadNeg=0.0;
		1434	flgstop=Cflgstop=0;
		1435	flagdeg=Cflagdeg=0;
		1436	flagstopdeg=Cflagstopdeg=0;
		1437	ekt=Cekt=0.0;
		1438	interup=Cinterup=0;
		1439	kcode=Ckcode=0;
		1440	kon=Ckon=0;
		1441	iCapNeg=CiCapNeg=0;
		1442	iNoFneg=CiNoFneg=0;
		1443	iNoFpos=CiNoFpos=0;
		1444	iCapPos=CiCapPos=0;
		1445	iDeg=CiDeg=0;
		1446	LP=CLP=0;
		1447	LN=CLN=0;
		1448	capSlope=CcapSlope=0.0;
		1449	capDispPos=CcapDispPos=0.0;
		1450	capDispNeg=CcapDispNeg=0.0;
		1451	elstk=Celstk=0.0;
		1452	fyieldPos=CfyieldPos=0.0;
		1453	fyieldNeg=CfyieldNeg=0.0;
		1454	alpha=Calpha=0.0;
		1455	ecaps=Cecaps=0.0;
		1456	ecapk=Cecapk=0.0;
		1457	ecapd=Cecapd=0.0;
		1458	cs=Ccs=0.0;
		1459	ck=Cck=0.0;
		1460	cd=Ccd=0.0;
		1461	dmax=Cdmax=0.0;
		1462	dmin=Cdmin=0.0;
		1463	Enrgtot=CEnrgtot=0.0;
		1464	Enrgc=CEnrgc=0.0;
		1465	fyPos=CfyPos=0.0;
		1466	fLimNeg=CfLimNeg=0.0;
		1467	fyNeg=CfyNeg=0.0;
		1468	ekP=CekP=0.0;
		1469	ekunload=Cekunload=0.0;
		1470	sp=Csp=0.0;
		1471	sn=Csn=0.0;
		1472	dP=CdP=0.0;
		1473	fP=CfP=0.0;
		1474	ek=Cek=0.0;
		1475	stif=Cstif=0.0;
		1476	dLimPos=CdLimPos=0.0;
		1477	dLimNeg=CdLimNeg=0.0;
		1478	vtot=Cvtot=0.0;
		1479	ftot=Cftot=0.0;
		1480	dtot=Cdtot=0.0;
		1481	dn=Cdn=0.0;
		1482	cpPos=CcpPos=0.0;
		1483	cpNeg=CcpNeg=0.0;
		1484	fLimPos=CfLimPos=0.0;
		1485	dlstPos=CdlstPos=0.0;
		1486	flstPos=CflstPos=0.0;
		1487	dlstNeg=CdlstNeg=0.0;
		1488	flstNeg=CflstNeg=0.0;
		1489	ekexcurs=Cekexcurs=0.0;
		1490	RSE=CRSE=0.0;
		1491	fPeakPos=CfPeakPos=0.0;
		1492	fPeakNeg=CfPeakNeg=0.0;
		1493	dCap1Pos=CdCap1Pos=0.0;
		1494	dCap2Pos=CdCap2Pos=0.0;
		1495	dCap1Neg=CdCap1Neg=0.0;
		1496	dCap2Neg=CdCap2Neg=0.0;
		1497	alphaNeg=CalphaNeg=0.0;
		1498	alphaPos=CalphaPos=0.0;
		1499	ekhardNeg=CekhardNeg=0.0;
		1500	ekhardPos=CekhardPos=0.0;
		1501	fCapRefPos=CfCapRefPos=0.0;
		1502	fCapRefNeg=CfCapRefNeg=0.0;
		1503	Enrgts=CEnrgts=0.0;
		1504	Enrgtk=CEnrgtk=0.0;
		1505	Enrgtd=CEnrgtd=0.0;
		1506	dyPos=CdyPos=0.0;
		1507	dyNeg=CdyNeg=0.0;
		1508	dyieldPos=CdyieldPos=0.0;
		1509	dyieldNeg=CdyieldNeg=0.0;
		1510	resSnHor=CresSnHor=0.0;
		1511	fmax=Cfmax=0.0;
		1512	fmin=Cfmin=0.0;
		1513	resSp=CresSp=0.0;
		1514	resSn=CresSn=0.0;
		1515	fCapPos=CfCapPos=0.0;
		1516	fCapNeg=CfCapNeg=0.0;
		1517	snHor=CsnHor=0.0;
		1518	spHor=CspHor=0.0;
		1519	resSpHor=CresSpHor=0.0;
		1520	snEnv=CsnEnv=0.0;
		1521	resSnEnv=CresSnEnv=0.0;
		1522	spEnv=CspEnv=0.0;
		1523	resSpEnv=CresSpEnv=0.0;
		1524	Resfac=CResfac=0.0;
		1525	capSlopeOrig=CcapSlopeOrig=0.0;
		1526	fracDispPos=CfracDispPos=0.0;
		1527	fracDispNeg=CfracDispNeg=0.0;
		1528	DPlus=CDPlus=0.0;
		1529	DNeg=CDNeg=0.0;
		1530	alphaN=CalphaN=0.0;
		1531	ecapa=Cecapa=0.0;
		1532	ca=Cca=0.0;
		1533	ResfacNeg=CResfacNeg=0.0;
		1534	myFcapping=CmyFcapping=0.0;
		1535	myFcappingNeg=CmyFcappingNeg=0.0;
		1536	capSlopeNeg=CcapSlopeNeg=0.0;
		1537	flagControlResponse=CflagControlResponse=0;
		1538	capSlopeOrigNeg=CcapSlopeOrigNeg=0.0;
		1539	Uprev=CUprev=0.0;
		1540	//then I initialize everything accordingly
		1541	elstk = Ke;
		1542	fyieldPos = My_pos;
		1543	fyieldNeg = My_neg;
		1544	alpha = As;
		1545	alphaN = AsNeg;
		1546	ecaps = LamdaS/(fyieldPos/(elstk));
		1547	ecapk = LamdaK/(fyieldPos/(elstk));
		1548	ecapa = LamdaA/(fyieldPos/(elstk));
		1549	ecapd = LamdaD/(fyieldPos/(elstk));
		1550	cs = Cs;
		1551	ck = Ck;
		1552	ca = Ca;
		1553	cd = Cd;
		1554	capDispPos = Thetap_pos+fyieldPos/elstk;
		1555	capDispNeg = -Thetap_neg+fyieldNeg/elstk;
		1556	Resfac = K;
		1557	ResfacNeg = KNeg;
		1558	myFcapping =-(fyieldPos+alphaelstkcapDispPos);
		1559	capSlope = myFcapping/(Thetapc_pos*elstk);
		1560	myFcappingNeg = -(abs(fyieldNeg)+alphaelstkabs(capDispNeg));
		1561	capSlopeNeg = myFcappingNeg/(Thetapc_neg*elstk);
		1562	fracDispPos = Thetau_pos;
		1563	fracDispNeg = -Thetau_neg;
		1564	DPlus = PDPlus;
		1565	DNeg = PDNeg;
		1566	stif = elstk;
		1567	ekP = elstk;
		1568	flagControlResponse = 0;
		1569	Tangent=Ke;
		1570	return 0;
		1571	}
		1572
		1573	UniaxialMaterial *
		1574	Bilin::getCopy(void)
		1575	{
		1576	Bilin *theCopy = new Bilin(this->getTag(),Ke,As,AsNeg,My_pos,My_neg,LamdaS,LamdaK,
		1577	LamdaA,LamdaD,Cs,Ck,Ca,Cd,Thetap_pos,Thetap_neg,Thetapc_pos,Thetapc_neg,K,KNeg,
		1578	Thetau_pos,Thetau_neg,PDPlus,PDNeg);
		1579
		1580	//Fixed Model parameters: need to change according to material properties
		1581	theCopy->U=U;
		1582	theCopy->CU=CU;
		1583	theCopy->Force=Force;
		1584	theCopy->CForce=CForce;
		1585	theCopy->Tangent=Tangent;
		1586	theCopy->CTangent=CTangent;
		1587	theCopy->dNewLoadPos=dNewLoadPos;
		1588	theCopy->CdNewLoadPos=CdNewLoadPos;
		1589	theCopy->dNewLoadNeg=dNewLoadNeg;
		1590	theCopy->CdNewLoadNeg=CdNewLoadNeg;
		1591	theCopy->flgstop=flgstop;
		1592	theCopy->Cflgstop=Cflgstop;
		1593	theCopy->flagdeg=flagdeg;
		1594	theCopy->Cflagdeg=Cflagdeg;
		1595	theCopy->flagstopdeg=flagstopdeg;
		1596	theCopy->Cflagstopdeg=Cflagstopdeg;
		1597	theCopy->ekt=ekt;
		1598	theCopy->Cekt=Cekt;
		1599	theCopy->interup=interup;
		1600	theCopy->Cinterup=Cinterup;
		1601	theCopy->kcode=kcode;
		1602	theCopy->Ckcode=Ckcode;
		1603	theCopy->kon=kon;
		1604	theCopy->Ckon=Ckon;
		1605	theCopy->iCapNeg=iCapNeg;
		1606	theCopy->CiCapNeg=CiCapNeg;
		1607	theCopy->iNoFneg=iNoFneg;
		1608	theCopy->CiNoFneg=CiNoFneg;
		1609	theCopy->iNoFpos=iNoFpos;
		1610	theCopy->CiNoFpos=CiNoFpos;
		1611	theCopy->iCapPos=iCapPos;
		1612	theCopy->CiCapPos=CiCapPos;
		1613	theCopy->iDeg=iDeg;
		1614	theCopy->CiDeg=CiDeg;
		1615	theCopy->LP=LP;
		1616	theCopy->CLP=CLP;
		1617	theCopy->LN=LN;
		1618	theCopy->CLN=CLN;
		1619	theCopy->capSlope=capSlope;
		1620	theCopy->CcapSlope=CcapSlope;
		1621	theCopy->capDispPos=capDispPos;
		1622	theCopy->CcapDispPos=CcapDispPos;
		1623	theCopy->capDispNeg=capDispNeg;
		1624	theCopy->CcapDispNeg=CcapDispNeg;
		1625	theCopy->elstk=elstk;
		1626	theCopy->Celstk=Celstk;
		1627	theCopy->fyieldPos=fyieldPos;
		1628	theCopy->CfyieldPos=CfyieldPos;
		1629	theCopy->fyieldNeg=fyieldNeg;
		1630	theCopy->CfyieldNeg=CfyieldNeg;
		1631	theCopy->alpha=alpha;
		1632	theCopy->Calpha=Calpha;
		1633	theCopy->ecaps=ecaps;
		1634	theCopy->Cecaps=Cecaps;
		1635	theCopy->ecapk=ecapk;
		1636	theCopy->Cecapk=Cecapk;
		1637	theCopy->ecapd=ecapd;
		1638	theCopy->Cecapd=Cecapd;
		1639	theCopy->cs=cs;
		1640	theCopy->Ccs=Ccs;
		1641	theCopy->ck=ck;
		1642	theCopy->Cck=Cck;
		1643	theCopy->cd=cd;
		1644	theCopy->Ccd=Ccd;
		1645	theCopy->dmax=dmax;
		1646	theCopy->Cdmax=Cdmax;
		1647	theCopy->dmin=dmin;
		1648	theCopy->Cdmin=Cdmin;
		1649	theCopy->Enrgtot=Enrgtot;
		1650	theCopy->CEnrgtot=CEnrgtot;
		1651	theCopy->Enrgc=Enrgc;
		1652	theCopy->CEnrgc=CEnrgc;
		1653	theCopy->fyPos=fyPos;
		1654	theCopy->CfyPos=CfyPos;
		1655	theCopy->fLimNeg=fLimNeg;
		1656	theCopy->CfLimNeg=CfLimNeg;
		1657	theCopy->fyNeg=fyNeg;
		1658	theCopy->CfyNeg=CfyNeg;
		1659	theCopy->ekP=ekP;
		1660	theCopy->CekP=CekP;
		1661	theCopy->ekunload=ekunload;
		1662	theCopy->Cekunload=Cekunload;
		1663	theCopy->sp=sp;
		1664	theCopy->Csp=Csp;
		1665	theCopy->sn=sn;
		1666	theCopy->Csn=Csn;
		1667	theCopy->dP=dP;
		1668	theCopy->CdP=CdP;
		1669	theCopy->fP=fP;
		1670	theCopy->CfP=CfP;
		1671	theCopy->ek=ek;
		1672	theCopy->Cek=Cek;
		1673	theCopy->stif=stif;
		1674	theCopy->Cstif=Cstif;
		1675	theCopy->dLimPos=dLimPos;
		1676	theCopy->CdLimPos=CdLimPos;
		1677	theCopy->dLimNeg=dLimNeg;
		1678	theCopy->CdLimNeg=CdLimNeg;
		1679	theCopy->vtot=vtot;
		1680	theCopy->Cvtot=Cvtot;
		1681	theCopy->ftot=ftot;
		1682	theCopy->Cftot=Cftot;
		1683	theCopy->dtot=dtot;
		1684	theCopy->Cdtot=Cdtot;
		1685	theCopy->dn=dn;
		1686	theCopy->Cdn=Cdn;
		1687	theCopy->cpPos=cpPos;
		1688	theCopy->CcpPos=CcpPos;
		1689	theCopy->cpNeg=cpNeg;
		1690	theCopy->CcpNeg=CcpNeg;
		1691	theCopy->fLimPos=fLimPos;
		1692	theCopy->CfLimPos=CfLimPos;
		1693	theCopy->dlstPos=dlstPos;
		1694	theCopy->CdlstPos=CdlstPos;
		1695	theCopy->flstPos=flstPos;
		1696	theCopy->CflstPos=CflstPos;
		1697	theCopy->dlstNeg=dlstNeg;
		1698	theCopy->CdlstNeg=CdlstNeg;
		1699	theCopy->flstNeg=flstNeg;
		1700	theCopy->CflstNeg=CflstNeg;
		1701	theCopy->ekexcurs=ekexcurs;
		1702	theCopy->Cekexcurs=Cekexcurs;
		1703	theCopy->RSE=RSE;
		1704	theCopy->CRSE=CRSE;
		1705	theCopy->fPeakPos=fPeakPos;
		1706	theCopy->CfPeakPos=CfPeakPos;
		1707	theCopy->fPeakNeg=fPeakNeg;
		1708	theCopy->CfPeakNeg=CfPeakNeg;
		1709	theCopy->dCap1Pos=dCap1Pos;
		1710	theCopy->CdCap1Pos=CdCap1Pos;
		1711	theCopy->dCap2Pos=dCap2Pos;
		1712	theCopy->CdCap2Pos=CdCap2Pos;
		1713	theCopy->dCap1Neg=dCap1Neg;
		1714	theCopy->CdCap1Neg=CdCap1Neg;
		1715	theCopy->dCap2Neg=dCap2Neg;
		1716	theCopy->CdCap2Neg=CdCap2Neg;
		1717	theCopy->alphaNeg=alphaNeg;
		1718	theCopy->CalphaNeg=CalphaNeg;
		1719	theCopy->alphaPos=alphaPos;
		1720	theCopy->CalphaPos=CalphaPos;
		1721	theCopy->ekhardNeg=ekhardNeg;
		1722	theCopy->CekhardNeg=CekhardNeg;
		1723	theCopy->ekhardPos=ekhardPos;
		1724	theCopy->CekhardPos=CekhardPos;
		1725	theCopy->fCapRefPos=fCapRefPos;
		1726	theCopy->CfCapRefPos=CfCapRefPos;
		1727	theCopy->fCapRefNeg=fCapRefNeg;
		1728	theCopy->CfCapRefNeg=CfCapRefNeg;
		1729	theCopy->Enrgts=Enrgts;
		1730	theCopy->CEnrgts=CEnrgts;
		1731	theCopy->Enrgtk=Enrgtk;
		1732	theCopy->CEnrgtk=CEnrgtk;
		1733	theCopy->Enrgtd=Enrgtd;
		1734	theCopy->CEnrgtd=CEnrgtd;
		1735	theCopy->dyPos=dyPos;
		1736	theCopy->CdyPos=CdyPos;
		1737	theCopy->dyNeg=dyNeg;
		1738	theCopy->CdyNeg=CdyNeg;
		1739	theCopy->dyieldPos=dyieldPos;
		1740	theCopy->CdyieldPos=CdyieldPos;
		1741	theCopy->dyieldNeg=dyieldNeg;
		1742	theCopy->CdyieldNeg=CdyieldNeg;
		1743	theCopy->resSnHor=resSnHor;
		1744	theCopy->CresSnHor=CresSnHor;
		1745	theCopy->fmax=fmax;
		1746	theCopy->Cfmax=Cfmax;
		1747	theCopy->fmin=fmin;
		1748	theCopy->Cfmin=Cfmin;
		1749	theCopy->resSp=resSp;
		1750	theCopy->CresSp=CresSp;
		1751	theCopy->resSn=resSn;
		1752	theCopy->CresSn=CresSn;
		1753	theCopy->fCapPos=fCapPos;
		1754	theCopy->CfCapPos=CfCapPos;
		1755	theCopy->fCapNeg=fCapNeg;
		1756	theCopy->CfCapNeg=CfCapNeg;
		1757	theCopy->snHor=snHor;
		1758	theCopy->CsnHor=CsnHor;
		1759	theCopy->spHor=spHor;
		1760	theCopy->CspHor=CspHor;
		1761	theCopy->resSpHor=resSpHor;
		1762	theCopy->CresSpHor=CresSpHor;
		1763	theCopy->snEnv=snEnv;
		1764	theCopy->CsnEnv=CsnEnv;
		1765	theCopy->resSnEnv=resSnEnv;
		1766	theCopy->CresSnEnv=CresSnEnv;
		1767	theCopy->spEnv=spEnv;
		1768	theCopy->CspEnv=CspEnv;
		1769	theCopy->resSpEnv=resSpEnv;
		1770	theCopy->CresSpEnv=CresSpEnv;
		1771	theCopy->Resfac=Resfac;
		1772	theCopy->CResfac=CResfac;
		1773	theCopy->capSlopeOrig=capSlopeOrig;
		1774	theCopy->CcapSlopeOrig=CcapSlopeOrig;
		1775	theCopy->fracDispPos=fracDispPos;
		1776	theCopy->CfracDispPos=CfracDispPos;
		1777	theCopy->fracDispNeg=fracDispNeg;
		1778	theCopy->CfracDispNeg=CfracDispNeg;
		1779	theCopy->DPlus=DPlus;
		1780	theCopy->CDPlus=CDPlus;
		1781	theCopy->DNeg=DNeg;
		1782	theCopy->CDNeg=CDNeg;
		1783	theCopy->alphaN=alphaN;
		1784	theCopy->CalphaN=CalphaN;
		1785	theCopy->ecapa=ecapa;
		1786	theCopy->Cecapa=Cecapa;
		1787	theCopy->ca=ca;
		1788	theCopy->Cca=Cca;
		1789	theCopy->ResfacNeg=ResfacNeg;
		1790	theCopy->CResfacNeg=CResfacNeg;
		1791	theCopy->myFcapping=myFcapping;
		1792	theCopy->CmyFcapping=CmyFcapping;
		1793	theCopy->myFcappingNeg=myFcappingNeg;
		1794	theCopy->CmyFcappingNeg=CmyFcappingNeg;
		1795	theCopy->capSlopeNeg=capSlopeNeg;
		1796	theCopy->CcapSlopeNeg=CcapSlopeNeg;
		1797	theCopy->flagControlResponse=flagControlResponse;
		1798	theCopy->CflagControlResponse=CflagControlResponse;
		1799	theCopy->capSlopeOrigNeg=capSlopeOrigNeg;
		1800	theCopy->CcapSlopeOrigNeg=CcapSlopeOrigNeg;
		1801	theCopy->CUprev=CUprev;
		1802	theCopy->Uprev=Uprev;
		1803
		1804	return theCopy;
		1805	}
		1806
		1807	int
		1808	Bilin::sendSelf(int cTag, Channel &theChannel)
		1809	{
		1810	int res = 0;
		1811
		1812	static Vector data(246); ///diorthose to megethos
		1813	data(0) = this->getTag();
		1814	data(1)=Ke;
		1815	data(2)=As;
		1816	data(3)=AsNeg;
		1817	data(4)=My_pos;
		1818	data(5)=My_neg;
		1819	data(6)=LamdaS;
		1820	data(7)=LamdaK;
		1821	data(8)=LamdaA;
		1822	data(9)=LamdaD;
		1823	data(10)=Cs;
		1824	data(11)=Ck;
		1825	data(12)=Ca;
		1826	data(13)=Cd;
		1827	data(14)=Thetap_pos;
		1828	data(15)=Thetap_neg;
		1829	data(16)=Thetapc_pos;
		1830	data(17)=Thetapc_neg;
		1831	data(18)=K;
		1832	data(19)=KNeg;
		1833	data(20)=Thetau_pos;
		1834	data(21)=Thetau_neg;
		1835	data(22)=PDPlus;
		1836	data(23)=PDNeg;
		1837	data(24)=U;
		1838	data(25)=CU;
		1839	data(26)=Force;
		1840	data(27)=CForce;
		1841	data(28)=Tangent;
		1842	data(29)=CTangent;
		1843	data(30)=dNewLoadPos;
		1844	data(31)=CdNewLoadPos;
		1845	data(32)=dNewLoadNeg;
		1846	data(33)=CdNewLoadNeg;
		1847	data(34)=flgstop;
		1848	data(35)=Cflgstop;
		1849	data(36)=flagdeg;
		1850	data(37)=Cflagdeg;
		1851	data(38)=flagstopdeg;
		1852	data(39)=Cflagstopdeg;
		1853	data(40)=ekt;
		1854	data(41)=Cekt;
		1855	data(42)=interup;
		1856	data(43)=Cinterup;
		1857	data(44)=kcode;
		1858	data(45)=Ckcode;
		1859	data(46)=kon;
		1860	data(47)=Ckon;
		1861	data(48)=iCapNeg;
		1862	data(49)=CiCapNeg;
		1863	data(50)=iNoFneg;
		1864	data(51)=CiNoFneg;
		1865	data(52)=iNoFpos;
		1866	data(53)=CiNoFpos;
		1867	data(54)=iCapPos;
		1868	data(55)=CiCapPos;
		1869	data(56)=iDeg;
		1870	data(57)=CiDeg;
		1871	data(58)=LP;
		1872	data(59)=CLP;
		1873	data(60)=LN;
		1874	data(61)=CLN;
		1875	data(62)=capSlope;
		1876	data(63)=CcapSlope;
		1877	data(64)=capDispPos;
		1878	data(65)=CcapDispPos;
		1879	data(66)=capDispNeg;
		1880	data(67)=CcapDispNeg;
		1881	data(68)=elstk;
		1882	data(69)=Celstk;
		1883	data(70)=fyieldPos;
		1884	data(71)=CfyieldPos;
		1885	data(72)=fyieldNeg;
		1886	data(73)=CfyieldNeg;
		1887	data(74)=alpha;
		1888	data(75)=Calpha;
		1889	data(76)=ecaps;
		1890	data(77)=Cecaps;
		1891	data(78)=ecapk;
		1892	data(79)=Cecapk;
		1893	data(80)=ecapd;
		1894	data(81)=Cecapd;
		1895	data(82)=cs;
		1896	data(83)=Ccs;
		1897	data(84)=ck;
		1898	data(85)=Cck;
		1899	data(86)=cd;
		1900	data(87)=Ccd;
		1901	data(88)=dmax;
		1902	data(89)=Cdmax;
		1903	data(90)=dmin;
		1904	data(91)=Cdmin;
		1905	data(92)=Enrgtot;
		1906	data(93)=CEnrgtot;
		1907	data(94)=Enrgc;
		1908	data(95)=CEnrgc;
		1909	data(96)=fyPos;
		1910	data(97)=CfyPos;
		1911	data(98)=fLimNeg;
		1912	data(99)=CfLimNeg;
		1913	data(100)=fyNeg;
		1914	data(101)=CfyNeg;
		1915	data(102)=ekP;
		1916	data(103)=CekP;
		1917	data(104)=ekunload;
		1918	data(105)=Cekunload;
		1919	data(106)=Csp;
		1920	data(107)=sn;
		1921	data(108)=Csn;
		1922	data(109)=dP;
		1923	data(110)=CdP;
		1924	data(111)=fP;
		1925	data(112)=CfP;
		1926	data(113)=ek;
		1927	data(114)=Cek;
		1928	data(115)=stif;
		1929	data(116)=Cstif;
		1930	data(117)=dLimPos;
		1931	data(118)=CdLimPos;
		1932	data(119)=dLimNeg;
		1933	data(120)=CdLimNeg;
		1934	data(121)=vtot;
		1935	data(122)=Cvtot;
		1936	data(123)=ftot;
		1937	data(124)=Cftot;
		1938	data(125)=dtot;
		1939	data(126)=Cdtot;
		1940	data(127)=dn;
		1941	data(128)=Cdn;
		1942	data(129)=cpPos;
		1943	data(130)=CcpPos;
		1944	data(131)=cpNeg;
		1945	data(132)=CcpNeg;
		1946	data(133)=fLimPos;
		1947	data(134)=CfLimPos;
		1948	data(135)=dlstPos;
		1949	data(136)=CdlstPos;
		1950	data(137)=flstPos;
		1951	data(138)=CflstPos;
		1952	data(139)=dlstNeg;
		1953	data(140)=CdlstNeg;
		1954	data(141)=flstNeg;
		1955	data(142)=CflstNeg;
		1956	data(143)=ekexcurs;
		1957	data(144)=Cekexcurs;
		1958	data(145)=RSE;
		1959	data(146)=CRSE;
		1960	data(147)=fPeakPos;
		1961	data(148)=CfPeakPos;
		1962	data(149)=fPeakNeg;
		1963	data(150)=CfPeakNeg;
		1964	data(151)=dCap1Pos;
		1965	data(152)=CdCap1Pos;
		1966	data(153)=dCap2Pos;
		1967	data(154)=CdCap2Pos;
		1968	data(155)=dCap1Neg;
		1969	data(156)=CdCap1Neg;
		1970	data(157)=dCap2Neg;
		1971	data(158)=CdCap2Neg;
		1972	data(159)=alphaNeg;
		1973	data(160)=CalphaNeg;
		1974	data(161)=alphaPos;
		1975	data(162)=CalphaPos;
		1976	data(163)=ekhardNeg;
		1977	data(164)=CekhardNeg;
		1978	data(165)=ekhardPos;
		1979	data(166)=CekhardPos;
		1980	data(167)=fCapRefPos;
		1981	data(168)=CfCapRefPos;
		1982	data(169)=fCapRefNeg;
		1983	data(170)=CfCapRefNeg;
		1984	data(171)=Enrgts;
		1985	data(172)=CEnrgts;
		1986	data(173)=Enrgtk;
		1987	data(174)=CEnrgtk;
		1988	data(175)=Enrgtd;
		1989	data(176)=CEnrgtd;
		1990	data(177)=dyPos;
		1991	data(178)=CdyPos;
		1992	data(179)=dyNeg;
		1993	data(180)=CdyNeg;
		1994	data(181)=dyieldPos;
		1995	data(182)=CdyieldPos;
		1996	data(183)=dyieldNeg;
		1997	data(184)=CdyieldNeg;
		1998	data(185)=resSnHor;
		1999	data(186)=CresSnHor;
		2000	data(187)=fmax;
		2001	data(188)=Cfmax;
		2002	data(189)=fmin;
		2003	data(190)=Cfmin;
		2004	data(191)=resSp;
		2005	data(192)=CresSp;
		2006	data(193)=resSn;
		2007	data(194)=CresSn;
		2008	data(195)=fCapPos;
		2009	data(196)=CfCapPos;
		2010	data(197)=fCapNeg;
		2011	data(198)=CfCapNeg;
		2012	data(199)=snHor;
		2013	data(200)=CsnHor;
		2014	data(201)=spHor;
		2015	data(202)=CspHor;
		2016	data(203)=resSpHor;
		2017	data(204)=CresSpHor;
		2018	data(205)=snEnv;
		2019	data(206)=CsnEnv;
		2020	data(207)=resSnEnv;
		2021	data(208)=CresSnEnv;
		2022	data(209)=spEnv;
		2023	data(210)=CspEnv;
		2024	data(211)=resSpEnv;
		2025	data(212)=CresSpEnv;
		2026	data(213)=Resfac;
		2027	data(214)=CResfac;
		2028	data(215)=capSlopeOrig;
		2029	data(216)=CcapSlopeOrig;
		2030	data(217)=fracDispPos;
		2031	data(218)=CfracDispPos;
		2032	data(219)=fracDispNeg;
		2033	data(220)=CfracDispNeg;
		2034	data(221)=DPlus;
		2035	data(222)=CDPlus;
		2036	data(223)=DNeg;
		2037	data(224)=CDNeg;
		2038	data(225)=alphaN;
		2039	data(226)=CalphaN;
		2040	data(227)=ecapa;
		2041	data(228)=Cecapa;
		2042	data(229)=ca;
		2043	data(230)=Cca;
		2044	data(231)=ResfacNeg;
		2045	data(232)=CResfacNeg;
		2046	data(233)=myFcapping;
		2047	data(234)=CmyFcapping;
		2048	data(235)=myFcappingNeg;
		2049	data(236)=CmyFcappingNeg;
		2050	data(237)=capSlopeNeg;
		2051	data(238)=CcapSlopeNeg;
		2052	data(239)=flagControlResponse;
		2053	data(240)=CflagControlResponse;
		2054	data(241)=sp;
		2055	data(242)=CcapSlopeOrigNeg;
		2056	data(243)=capSlopeOrigNeg;
		2057	data(244)=CUprev;
		2058	data(245)=Uprev;
		2059
		2060	res = theChannel.sendVector(this->getDbTag(), cTag, data);
		2061	if (res < 0)
		2062	opserr << "Bilin::sendSelf() - failed to send data\n";
		2063
		2064	return res;
		2065	}
		2066
		2067	int
		2068	Bilin::recvSelf(int cTag, Channel &theChannel,
		2069	FEM_ObjectBroker &theBroker)
		2070	{
		2071	int res = 0;
		2072	static Vector data(246);
		2073	res = theChannel.recvVector(this->getDbTag(), cTag, data);
		2074
		2075	if (res < 0) {
		2076	opserr << "Bilin::recvSelf() - failed to receive data\n";
		2077	this->setTag(0);
		2078	}
		2079	else {
		2080	this->setTag((int)data(0));
		2081	Ke=data(1);
		2082	As=data(2);
		2083	AsNeg=data(3);
		2084	My_pos=data(4);
		2085	My_neg=data(5);
		2086	LamdaS=data(6);
		2087	LamdaK=data(7);
		2088	LamdaA=data(8);
		2089	LamdaD=data(9);
		2090	Cs=data(10);
		2091	Ck=data(11);
		2092	Ca=data(12);
		2093	Cd=data(13);
		2094	Thetap_pos=data(14);
		2095	Thetap_neg=data(15);
		2096	Thetapc_pos=data(16);
		2097	Thetapc_neg=data(17);
		2098	K=data(18);
		2099	KNeg=data(19);
		2100	Thetau_pos=data(20);
		2101	Thetau_neg=data(21);
		2102	PDPlus=data(22);
		2103	PDNeg=data(23);
		2104	U=data(24);
		2105	CU=data(25);
		2106	Force=data(26);
		2107	CForce=data(27);
		2108	Tangent=data(28);
		2109	CTangent=data(29);
		2110	dNewLoadPos=data(30);
		2111	CdNewLoadPos=data(31);
		2112	dNewLoadNeg=data(32);
		2113	CdNewLoadNeg=data(33);
		2114	flgstop=data(34);
		2115	Cflgstop=data(35);
		2116	flagdeg=data(36);
		2117	Cflagdeg=data(37);
		2118	flagstopdeg=data(38);
		2119	Cflagstopdeg=data(39);
		2120	ekt=data(40);
		2121	Cekt=data(41);
		2122	interup=data(42);
		2123	Cinterup=data(43);
		2124	kcode=data(44);
		2125	Ckcode=data(45);
		2126	kon=data(46);
		2127	Ckon=data(47);
		2128	iCapNeg=data(48);
		2129	CiCapNeg=data(49);
		2130	iNoFneg=data(50);
		2131	CiNoFneg=data(51);
		2132	iNoFpos=data(52);
		2133	CiNoFpos=data(53);
		2134	iCapPos=data(54);
		2135	CiCapPos=data(55);
		2136	iDeg=data(56);
		2137	CiDeg=data(57);
		2138	LP=data(58);
		2139	CLP=data(59);
		2140	LN=data(60);
		2141	CLN=data(61);
		2142	capSlope=data(62);
		2143	CcapSlope=data(63);
		2144	capDispPos=data(64);
		2145	CcapDispPos=data(65);
		2146	capDispNeg=data(66);
		2147	CcapDispNeg=data(67);
		2148	elstk=data(68);
		2149	Celstk=data(69);
		2150	fyieldPos=data(70);
		2151	CfyieldPos=data(71);
		2152	fyieldNeg=data(72);
		2153	CfyieldNeg=data(73);
		2154	alpha=data(74);
		2155	Calpha=data(75);
		2156	ecaps=data(76);
		2157	Cecaps=data(77);
		2158	ecapk=data(78);
		2159	Cecapk=data(79);
		2160	ecapd=data(80);
		2161	Cecapd=data(81);
		2162	cs=data(82);
		2163	Ccs=data(83);
		2164	ck=data(84);
		2165	Cck=data(85);
		2166	cd=data(86);
		2167	Ccd=data(87);
		2168	dmax=data(88);
		2169	Cdmax=data(89);
		2170	dmin=data(90);
		2171	Cdmin=data(91);
		2172	Enrgtot=data(92);
		2173	CEnrgtot=data(93);
		2174	Enrgc=data(94);
		2175	CEnrgc=data(95);
		2176	fyPos=data(96);
		2177	CfyPos=data(97);
		2178	fLimNeg=data(98);
		2179	CfLimNeg=data(99);
		2180	fyNeg=data(100);
		2181	CfyNeg=data(101);
		2182	ekP=data(102);
		2183	CekP=data(103);
		2184	ekunload=data(104);
		2185	Cekunload=data(105);
		2186	Csp=data(106);
		2187	sn=data(107);
		2188	Csn=data(108);
		2189	dP=data(109);
		2190	CdP=data(110);
		2191	fP=data(111);
		2192	CfP=data(112);
		2193	ek=data(113);
		2194	Cek=data(114);
		2195	stif=data(115);
		2196	Cstif=data(116);
		2197	dLimPos=data(117);
		2198	CdLimPos=data(118);
		2199	dLimNeg=data(119);
		2200	CdLimNeg=data(120);
		2201	vtot=data(121);
		2202	Cvtot=data(122);
		2203	ftot=data(123);
		2204	Cftot=data(124);
		2205	dtot=data(125);
		2206	Cdtot=data(126);
		2207	dn=data(127);
		2208	Cdn=data(128);
		2209	cpPos=data(129);
		2210	CcpPos=data(130);
		2211	cpNeg=data(131);
		2212	CcpNeg=data(132);
		2213	fLimPos=data(133);
		2214	CfLimPos=data(134);
		2215	dlstPos=data(135);
		2216	CdlstPos=data(136);
		2217	flstPos=data(137);
		2218	CflstPos=data(138);
		2219	dlstNeg=data(139);
		2220	CdlstNeg=data(140);
		2221	flstNeg=data(141);
		2222	CflstNeg=data(142);
		2223	ekexcurs=data(143);
		2224	Cekexcurs=data(144);
		2225	RSE=data(145);
		2226	CRSE=data(146);
		2227	fPeakPos=data(147);
		2228	CfPeakPos=data(148);
		2229	fPeakNeg=data(149);
		2230	CfPeakNeg=data(150);
		2231	dCap1Pos=data(151);
		2232	CdCap1Pos=data(152);
		2233	dCap2Pos=data(153);
		2234	CdCap2Pos=data(154);
		2235	dCap1Neg=data(155);
		2236	CdCap1Neg=data(156);
		2237	dCap2Neg=data(157);
		2238	CdCap2Neg=data(158);
		2239	alphaNeg=data(159);
		2240	CalphaNeg=data(160);
		2241	alphaPos=data(161);
		2242	CalphaPos=data(162);
		2243	ekhardNeg=data(163);
		2244	CekhardNeg=data(164);
		2245	ekhardPos=data(165);
		2246	CekhardPos=data(166);
		2247	fCapRefPos=data(167);
		2248	CfCapRefPos=data(168);
		2249	fCapRefNeg=data(169);
		2250	CfCapRefNeg=data(170);
		2251	Enrgts=data(171);
		2252	CEnrgts=data(172);
		2253	Enrgtk=data(173);
		2254	CEnrgtk=data(174);
		2255	Enrgtd=data(175);
		2256	CEnrgtd=data(176);
		2257	dyPos=data(177);
		2258	CdyPos=data(178);
		2259	dyNeg=data(179);
		2260	CdyNeg=data(180);
		2261	dyieldPos=data(181);
		2262	CdyieldPos=data(182);
		2263	dyieldNeg=data(183);
		2264	CdyieldNeg=data(184);
		2265	resSnHor=data(185);
		2266	CresSnHor=data(186);
		2267	fmax=data(187);
		2268	Cfmax=data(188);
		2269	fmin=data(189);
		2270	Cfmin=data(190);
		2271	resSp=data(191);
		2272	CresSp=data(192);
		2273	resSn=data(193);
		2274	CresSn=data(194);
		2275	fCapPos=data(195);
		2276	CfCapPos=data(196);
		2277	fCapNeg=data(197);
		2278	CfCapNeg=data(198);
		2279	snHor=data(199);
		2280	CsnHor=data(200);
		2281	spHor=data(201);
		2282	CspHor=data(202);
		2283	resSpHor=data(203);
		2284	CresSpHor=data(204);
		2285	snEnv=data(205);
		2286	CsnEnv=data(206);
		2287	resSnEnv=data(207);
		2288	CresSnEnv=data(208);
		2289	spEnv=data(209);
		2290	CspEnv=data(210);
		2291	resSpEnv=data(211);
		2292	CresSpEnv=data(212);
		2293	Resfac=data(213);
		2294	CResfac=data(214);
		2295	capSlopeOrig=data(215);
		2296	CcapSlopeOrig=data(216);
		2297	fracDispPos=data(217);
		2298	CfracDispPos=data(218);
		2299	fracDispNeg=data(219);
		2300	CfracDispNeg=data(220);
		2301	DPlus=data(221);
		2302	CDPlus=data(222);
		2303	DNeg=data(223);
		2304	CDNeg=data(224);
		2305	alphaN=data(225);
		2306	CalphaN=data(226);
		2307	ecapa=data(227);
		2308	Cecapa=data(228);
		2309	ca=data(229);
		2310	Cca=data(230);
		2311	ResfacNeg=data(231);
		2312	CResfacNeg=data(232);
		2313	myFcapping=data(233);
		2314	CmyFcapping=data(234);
		2315	myFcappingNeg=data(235);
		2316	CmyFcappingNeg=data(236);
		2317	capSlopeNeg=data(237);
		2318	CcapSlopeNeg=data(238);
		2319	flagControlResponse=data(239);
		2320	CflagControlResponse=data(240);
		2321	sp=data(241);
		2322	CcapSlopeOrigNeg=data(242);
		2323	capSlopeOrigNeg=data(243);
		2324	CUprev=data(244);
		2325	Uprev=data(245);
		2326	}
		2327
		2328	return res;
		2329	}
		2330
		2331	void
		2332	Bilin::Print(OPS_Stream &s, int flag)
		2333	{
		2334	//s << "Bilin tag: " << this->getTag() << endln;
		2335	//s << " G: " << G << endln;
		2336	//s << " t: " << t << endln;
		2337	//s << " A: " << A << endln;
		2338	}
		2339
		2340	//my functions
		2341	void
		2342	Bilin::interPoint(double &xInt, double &yInt, double x1, double y1, double m1, double x2, double y2, double m2)
		2343	{
		2344	xInt = (-m2x2+y2+m1x1-y1) / (m1-m2);
		2345	yInt = m1xInt-m1x1+y1;
		2346	}
		2347	void Bilin::snCalc(void)
		2348	{
		2349	// c Each time that the hysteretic loop changes from negative to positive
		2350	// c delta displacement, this subroutine calculates the point where ends
		2351	// c ekunload and starts whether the positive hardening curve or the
		2352	// c positive cap curve. In cases where "cpPos<fyPos" this point could
		2353	// c not be reached, because the unloading stiffness could intersect the
		2354	// c positive cap curve. However, this change is reflected in the main
		2355	// c program.
		2356	// c
		2357	// c Output Variables: sn,resSn,snHard,resSnHard
		2358	// c Input Variables: dP,fP,ekunload,alphaPos,dyPos,fyPos,cpPos,fCapPos
		2359	// c capStiff,fCapRefPos
		2360
		2361	double Resid = Resfac*fyPos;
		2362	double dresid = cpPos+(Resid-fCapPos)/(capSlope*elstk);
		2363	double ekresid = 1.0e-10;
		2364	dyPos = fyPos/elstk;
		2365	double snHard,resSnHard,snLim,resSnLim,snResid,resSnResid;
		2366	if (dyPos<cpPos) {
		2367	// call interPoint(snHard,resSnHard,dyPos,fyPos,elstk*alphaPos,dP,fP,ekunload)
		2368
		2369	interPoint(snHard,resSnHard,dyPos,fyPos,elstk*alphaPos,dP,fP,ekunload);
		2370	} else {
		2371	// call interPoint(snHard,resSnHard,cpPos,fCapPos,elstk*alphaPos,dP,fP,ekunload)
		2372
		2373	interPoint(snHard,resSnHard,cpPos,fCapPos,elstk*alphaPos,dP,fP,ekunload);
		2374	}
		2375
		2376	// call interPoint(snCap,resSnCap,0.d0,fCapRefPos,capSlope*elstk,dP,fP,ekunload)
		2377	double snCap,resSnCap;
		2378	interPoint(snCap,resSnCap,0.0,fCapRefPos,capSlope*elstk,dP,fP,ekunload);
		2379
		2380	//sn = min(snHard,snCap);
		2381	if (snHard<snCap)
		2382	{
		2383	sn=snHard;
		2384	}
		2385	else
		2386	{
		2387	sn=snCap;
		2388	}
		2389	//resSn = min(resSnHard,resSnCap);
		2390	if (resSnHard<resSnCap)
		2391	{
		2392	resSn=resSnHard;
		2393	}
		2394	else
		2395	{
		2396	resSn=resSnCap;
		2397	}
		2398	snEnv = sn;
		2399	resSnEnv = resSn;
		2400	// call interPoint(temp_1_farzin,temp,0.d0,fCapRefPos,capSlope*elstk,0.d0,Resid,0.d0)
		2401	//don't need to call that!!!!!!!!!!!!!!!
		2402	if((LP==1)&&(fLimPos==0.0)) {
		2403	//call interPoint(snLim,resSnLim,dLimPos,fLimPos,0.d0,dP,fP,ekunload)
		2404	interPoint(snLim,resSnLim,dLimPos,fLimPos,0.0,dP,fP,ekunload);
		2405	if (snLim<sn){
		2406	sn=snLim;
		2407	resSn=resSnLim;
		2408	}
		2409
		2410	//call interPoint(snHor,resSnHor,dLimPos,fLimPos,0.d0,dyPos,fyPos,elstk*alphaPos)
		2411	interPoint(snHor,resSnHor,dLimPos,fLimPos,0.0,dyPos,fyPos,elstk*alphaPos);
		2412	}
		2413
		2414	if (sn>dresid) {
		2415	// call interPoint(snResid,resSnResid,dresid,Resid,ekresid,dP,fP,ekunload)
		2416	interPoint(snResid,resSnResid,dresid,Resid,ekresid,dP,fP,ekunload);
		2417	sn = snResid;
		2418	resSn = resSnResid;
		2419	}
		2420	}
		2421
		2422
		2423	void
		2424	Bilin::envelPosCap2(double fy,double alphaPos,double alphaCap,double cpDsp,double& d,
		2425	double& f,double& ek,double elstk,double fyieldPos,double Resfac)
		2426	{
		2427
		2428
		2429	double fracDispPosV=0.20;
		2430	double dy = fy/elstk;
		2431
		2432	double Res,rcap,dres;
		2433	if(dy<=cpDsp) {
		2434	Res = Resfac*fyieldPos;
		2435	rcap = fy+alphaPoselstk(cpDsp-dy);
		2436	dres = cpDsp+(Res-rcap)/(alphaCap*elstk);
		2437
		2438	if (d<0.0){
		2439	f = 0.0;
		2440	ek = 1.0e-7;
		2441	}else if (d<=dy) {
		2442	ek = elstk;
		2443	f = ek*d;
		2444	} else if(d<=cpDsp) {
		2445	ek = elstk*alphaPos;
		2446	f = fy+ek*(d-dy);
		2447	} else if(d<=dres) {
		2448	ek = alphaCap*elstk;
		2449	f = rcap+ek*(d-cpDsp);
		2450	} else {
		2451	ek = 1.0e-7;
		2452	f = Res+d*ek;
		2453	}
		2454	// c added by Dimitrios to account for fracture
		2455	if(d>=fracDispPos) {
		2456	ek = 1.0e-7;
		2457	f = 1.0e-10;
		2458	d=fracDispPosV;
		2459	flgstop=1;
		2460	}
		2461	} else if(dy>cpDsp) {
		2462
		2463	rcap = elstk*cpDsp;
		2464	Res = Resfac*rcap;
		2465	dres = cpDsp+(Res-rcap)/(alphaCap*elstk);
		2466
		2467	if (d<0.0) {
		2468	f = 0.0;
		2469	ek = 1.0e-7;
		2470	} else if(d<=cpDsp) {
		2471	ek = elstk;
		2472	f = ek*d;
		2473	} else if(d<=dres) {
		2474	ek = alphaCap*elstk;
		2475	f = rcap+ek*(d-cpDsp);
		2476	} else {
		2477	ek = 1.0e-7;
		2478	f = Res+d*ek;
		2479	}
		2480	// c added by Dimitrios to account for fracture
		2481	if(d>=fracDispPos) {
		2482	ek = 1.0e-7;
		2483	f = 1.0e-10;
		2484	d=fracDispPosV;
		2485	flgstop=1;
		2486	}
		2487
		2488	} else {
		2489	}
		2490	}
		2491	double
		2492	Bilin::boundPos(void)
		2493	{
		2494
		2495	double Resid=0.0;
		2496	double dBoundPos;
		2497	dyNeg = fyNeg/elstk;
		2498	double dresid = cpPos+(Resid-fCapPos)/(capSlope*elstk);
		2499	double ekresid = 1.0e-10;
		2500
		2501	// call interPoint(d1,f1,E
		2502	double d1,f1;
		2503	interPoint(d1,f1,dyNeg,fyNeg,elstkalphaNeg,0.0,fCapRefPos,capSlopeelstk);
		2504
		2505	// call interPoint(d2,f2,dyNeg,fyNeg,elstk*alphaNeg,dresid,resid,ekresid)
		2506	double d2,f2;
		2507	interPoint(d2,f2,dyNeg,fyNeg,elstk*alphaNeg,dresid,Resid,ekresid);
		2508	//dBoundPos = max(d1,d2);
		2509	if (d1>d2)
		2510	{
		2511	dBoundPos=d1;
		2512	}
		2513	else
		2514	{
		2515	dBoundPos=d2;
		2516	}
		2517	return(dBoundPos);
		2518	}
		2519	void
		2520	Bilin::envHitsZero(double& f)
		2521	{
		2522	if (fP>0) {
		2523	if ((f*fP)<0) {
		2524	f = 0;
		2525	ek = 1.0e-7;
		2526	iNoFpos = 1;
		2527	flagControlResponse = 1;
		2528	}
		2529	}else if (fP<0) {
		2530	if ((f*fP)<0) {
		2531	f = 0;
		2532	ek = 1.0e-7;
		2533	iNoFneg = 1;
		2534	flagControlResponse = 1;
		2535	}
		2536	} else {
		2537	}
		2538	}
		2539	void
		2540	Bilin::envelNegCap2(double fy,double alphaNeg,double alphaCap,double cpDsp,double& d,double& f,double& ek,
		2541	double elstk,double fyieldNeg,double Resfac)
		2542	{
		2543
		2544	double fracDispNegV = -0.20;
		2545	double dy = fy/elstk;
		2546	double Res,rcap,dres;
		2547	if(dy>=cpDsp){
		2548
		2549	Res = Resfac*fyieldNeg;
		2550	rcap = fy+alphaNegelstk(cpDsp-dy);
		2551	dres = cpDsp+(Res-rcap)/(alphaCap*elstk);
		2552
		2553	if (d>0.0){
		2554	f = 0.0;
		2555	ek = 1.0e-7;
		2556	} else if (d>=dy) {
		2557	ek = elstk;
		2558	f = ek*d;
		2559	} else if (d>=cpDsp) {
		2560	ek = elstk*alphaNeg;
		2561	f = fy+ek*(d-dy);
		2562	} else if (d>=dres) {
		2563	ek = elstk*alphaCap;
		2564	f = rcap+ek*(d-cpDsp);
		2565	} else {
		2566	ek = 1.0e-7;
		2567	f = Res+ek*d;
		2568	}
		2569	//% c added by Dimitrios to account for fracture
		2570	if(d<=fracDispNegV) {
		2571	ek = 1.0e-7;
		2572	f = 1.0e-10;
		2573	d=fracDispNeg;
		2574	flgstop=1;
		2575	}
		2576
		2577	} else if(dy<cpDsp) {
		2578
		2579	rcap = elstk*cpDsp;
		2580	Res = Resfac*rcap;
		2581	dres = cpDsp+(Res-rcap)/(alphaCap*elstk);
		2582
		2583	if (d>0.0) {
		2584	f = 0.0;
		2585	ek = 1.0e-7;
		2586	} else if (d>=cpDsp) {
		2587	ek = elstk;
		2588	f = ek*d;
		2589	} else if (d>=dres) {
		2590	ek = elstk*alphaCap;
		2591	f = rcap+ek*(d-cpDsp);
		2592	} else {
		2593	ek = 1.0e-7;
		2594	f = Res+ek*d;
		2595	}
		2596	// c added by Dimitrios to account for fracture
		2597	if(d<=fracDispNegV) {
		2598	ek = 1.0e-7;
		2599	f = 1.0e-10;
		2600	d=fracDispNeg;
		2601	flgstop=1;
		2602	}
		2603
		2604	} else {
		2605	}
		2606	}
		2607	void
		2608	Bilin::spCalc(void)
		2609	{
		2610	double Resid = ResfacNeg*fyNeg;
		2611	dyNeg = fyNeg/elstk;
		2612	double dresid = cpNeg+(Resid-fCapNeg)/(capSlopeNeg*elstk);
		2613	double ekresid = 1.0e-10;
		2614	double spHard,resSpHard,spCap,resSpCap,spLim,resSpLim,spResid,resSpResid;
		2615	if (dyNeg>cpNeg) {
		2616	// call interPoint(spHard,resSpHard,dyNeg,fyNeg,elstk*alphaNeg,dP,fP,ekunload)
		2617	interPoint(spHard,resSpHard,dyNeg,fyNeg,elstk*alphaNeg,dP,fP,ekunload);
		2618	} else {
		2619	// call interPoint(spHard,resSpHard,cpNeg,fCapNeg,elstk*alphaNeg,dP,fP,ekunload)
		2620	interPoint(spHard,resSpHard,cpNeg,fCapNeg,elstk*alphaNeg,dP,fP,ekunload);
		2621	}
		2622
		2623	// call interpoint(spCap,resSpCap,0.d0,fCapRefNeg,capSlope*elstk,dP,fP,ekunload)
		2624	interPoint(spCap,resSpCap,0.0,fCapRefNeg,capSlopeNeg*elstk,dP,fP,ekunload);
		2625	//sp = max(spHard,spCap);
		2626	if(spHard>spCap)
		2627	{
		2628	sp=spHard;
		2629	}
		2630	else
		2631	{
		2632	sp=spCap;
		2633	}
		2634	//resSp = max(resSpHard,resSpCap);
		2635	if(resSpHard>resSpCap)
		2636	{
		2637	resSp=resSpHard;
		2638	}
		2639	else
		2640	{
		2641	resSp=resSpCap;
		2642	}
		2643	spEnv = sp;
		2644	resSpEnv = resSp;
		2645
		2646	if((LN==1)&&(fLimNeg==0.0)) {
		2647	// call interPoint(spLim,resSpLim,dLimNeg,fLimNeg,0.d0,dP,fP,ekunload)
		2648	interPoint(spLim,resSpLim,dLimNeg,fLimNeg,0.0,dP,fP,ekunload);
		2649	if (spLim>sp) {
		2650	sp=spLim;
		2651	resSp=resSpLim;
		2652	}
		2653
		2654	// call interPoint(spHor,resSpHor,dLimNeg,fLimNeg,0.d0,dyNeg,fyNeg,elstk*alphaNeg)
		2655	interPoint(spHor,resSpHor,dLimNeg,fLimNeg,0.0,dyNeg,fyNeg,elstk*alphaNeg);
		2656	}
		2657
		2658	if (sp<dresid) {
		2659	// call interPoint(spResid,resSpResid,dresid,Resid,ekresid,dP,fP,ekunload)
		2660	interPoint(spResid,resSpResid,dresid,Resid,ekresid,dP,fP,ekunload);
		2661	sp = spResid;
		2662	resSp = resSpResid;
		2663	}
		2664	}
		2665	double
		2666	Bilin::boundNeg(void)
		2667	{
		2668
		2669	double Resid = 0;
		2670	double dBoundNeg;
		2671	dyPos = fyPos/elstk;
		2672	double dresid = cpNeg+(Resid-fCapNeg)/(capSlopeNeg*elstk);
		2673	double ekresid = 1.0e-10;
		2674	double d1,f1,d2,f2;
		2675	// call interPoint(d1,f1,dyPos,fyPos,elstkalphaPos,0.d0,fCapRefNeg,capSlopeelstk)
		2676	interPoint(d1,f1,dyPos,fyPos,elstkalphaPos,0.0,fCapRefNeg,capSlopeNegelstk);
		2677	//%
		2678	// call interPoint(d2,f2,dyPos,fyPos,elstk*alphaPos,dresid,resid,ekresid)
		2679	interPoint(d2,f2,dyPos,fyPos,elstk*alphaPos,dresid,Resid,ekresid);
		2680	//dBoundNeg = min(d1,d2);
		2681	if (d1<d2)
		2682	{
		2683	dBoundNeg=d1;
		2684	}
		2685	else
		2686	{
		2687	dBoundNeg=d2;
		2688	}
		2689	return (dBoundNeg);
		2690	}

Subversion Repositories OpenSees

(root)/trunk/SRC/material/uniaxial/snap/Bilin.cpp - Rev 4147