Subversion Repositories OpenSees

Compare Revisions

Ignore whitespace Rev 6564 → Rev 6565

/trunk/SRC/element/Element.h
118,13 → 118,14
Matrix **previousK;
int numPreviousK;
 
private:
int index, nodeIndex;
 
static Matrix ** theMatrices;
static Vector ** theVectors1;
static Vector ** theVectors2;
static int numMatrices;
 
private:
};
 
 
/trunk/SRC/element/TclElementCommands.cpp
124,6 → 124,8
extern void *OPS_ShellMITC9(void);
extern void *OPS_ShellDKGQ(void); //Added by Lisha Wang, Xinzheng Lu, Linlin Xie, Song Cen & Quan Gu
extern void *OPS_ShellNLDKGQ(void); //Added by Lisha Wang, Xinzheng Lu, Linlin Xie, Song Cen & Quan Gu
extern void *OPS_ShellDKGT(void); //Added by Shuhao Zhang and Xinzheng Lu
extern void *OPS_ShellNLDKGT(void); //Added by Shuhao Zhang and Xinzheng Lu
extern void *OPS_Quad4FiberOverlay(void);
extern void *OPS_Brick8FiberOverlay(void);
extern void *OPS_QuadBeamEmbedContact(void);
709,8 → 711,28
opserr << "TclElementCommand -- unable to create element of type : " << argv[1] << endln;
return TCL_ERROR;
}
 
} else if ((strcmp(argv[1],"shellDKGT") == 0) || (strcmp(argv[1],"ShellDKGT") == 0)) {
void *theEle = OPS_ShellDKGT();
if (theEle != 0)
theElement = (Element *)theEle;
else {
opserr << "TclElementCommand -- unable to create element of type : " << argv[1] << endln;
return TCL_ERROR;
}
 
} else if ((strcmp(argv[1],"shellNLDKGT"G) == 0) || (strcmp(argv[1],"ShellNLDKGT") == 0)) {
void *theEle = OPS_ShellNLDKGT();
if (theEle != 0)
theElement = (Element *)theEle;
else {
opserr << "TclElementCommand -- unable to create element of type : " << argv[1] << endln;
return TCL_ERROR;
}
} else if ((strcmp(argv[1],"CoupledZeroLength") == 0) || (strcmp(argv[1],"ZeroLengthCoupled") == 0)) {
void *theEle = OPS_CoupledZeroLength();
/trunk/SRC/element/dispBeamColumn/DispBeamColumn3d.cpp
1253,8 → 1253,10
 
for (int i = 0; i < numSections; i++) {
opserr << "Section Type: " << theSections[i]->getClassTag() << endln;
theSections[i]->Print(s,flag);
// theSections[i]->Print(s,flag);
}
if (rho != 0)
opserr << "Mass: \n" << this->getMass();
}
 
 
/trunk/SRC/element/shell/ShellDKGT.h
196,4 → 196,4
// vector for applying loads
Vector *load;
Matrix *Ki;
} ;
} ;
/trunk/SRC/element/shell/ShellMITC4.cpp
67,7 → 67,7
int numArgs = OPS_GetNumRemainingInputArgs();
if (numArgs < 6) {
opserr << "Want: element ShellMITC4 $tag $iNode $jNoe $kNode $lNode $secTag";
opserr << "Want: element ShellMITC4 $tag $iNode $jNoe $kNode $lNode $secTag<-updateBasis>";
return 0;
}
77,7 → 77,14
opserr << "WARNING invalid integer tag: element ShellMITC4 \n";
return 0;
}
bool updateBasis = false;
 
if (numArgs == 7) {
const char* type = OPS_GetString();
if(strcmp(type,"-updateBasis") == 0)
updateBasis = true;
}
 
SectionForceDeformation *theSection = OPS_getSectionForceDeformation(iData[5]);
 
if (theSection == 0) {
86,7 → 93,7
}
theElement = new ShellMITC4(iData[0], iData[1], iData[2], iData[3],
iData[4], *theSection);
iData[4], *theSection, updateBasis);
 
return theElement;
}
110,7 → 117,7
//null constructor
ShellMITC4::ShellMITC4( ) :
Element( 0, ELE_TAG_ShellMITC4 ),
connectedExternalNodes(4), load(0), Ki(0)
connectedExternalNodes(4), load(0), Ki(0), doUpdateBasis(false)
{
for (int i = 0 ; i < 4; i++ )
materialPointers[i] = 0;
137,11 → 144,12
ShellMITC4::ShellMITC4( int tag,
int node1,
int node2,
int node3,
int node3,
int node4,
SectionForceDeformation &theMaterial ) :
SectionForceDeformation &theMaterial,
bool UpdateBasis) :
Element( tag, ELE_TAG_ShellMITC4 ),
connectedExternalNodes(4), load(0), Ki(0)
connectedExternalNodes(4), load(0), Ki(0), doUpdateBasis(UpdateBasis)
{
int i;
 
1126,7 → 1134,8
resid.Zero( ) ;
//start Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
updateBasis( );
if (doUpdateBasis == true)
updateBasis( );
//end Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
 
double dx34 = xl[0][2]-xl[0][3];
1891,7 → 1900,7
// Now quad sends the ids of its materials
int matDbTag;
static ID idData(13);
static ID idData(14);
int i;
for (i = 0; i < 4; i++) {
1912,6 → 1921,10
idData(10) = connectedExternalNodes(1);
idData(11) = connectedExternalNodes(2);
idData(12) = connectedExternalNodes(3);
if (doUpdateBasis == true)
idData(13) = 0;
else
idData(13) = 1;
 
res += theChannel.sendID(dataTag, commitTag, idData);
if (res < 0) {
1952,7 → 1965,7
int dataTag = this->getDbTag();
 
static ID idData(13);
static ID idData(14);
// Quad now receives the tags of its four external nodes
res += theChannel.recvID(dataTag, commitTag, idData);
if (res < 0) {
1965,6 → 1978,10
connectedExternalNodes(1) = idData(10);
connectedExternalNodes(2) = idData(11);
connectedExternalNodes(3) = idData(12);
if (idData(13) == 0)
doUpdateBasis = true;
else
doUpdateBasis = false;
 
static Vector vectData(5);
res += theChannel.recvVector(dataTag, commitTag, vectData);
/trunk/SRC/element/shell/ShellMITC4.h
57,7 → 57,8
int node2,
int node3,
int node4,
SectionForceDeformation &theMaterial ) ;
SectionForceDeformation &theMaterial,
bool updateBasis=false) ;
//destructor
virtual ~ShellMITC4( ) ;
150,9 → 151,10
 
//compute local coordinates and basis
void computeBasis( ) ;
//start Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
//start Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
bool doUpdateBasis;
void updateBasis( ) ;
//end Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
//end Yuli Huang (yulihuang@gmail.com) & Xinzheng Lu (luxz@tsinghua.edu.cn)
//inertia terms
void formInertiaTerms( int tangFlag ) ;
190,8 → 192,4
 
 
 
 
 
 
 
#endif
/trunk/SRC/element/componentElement/ComponentElement2d.cpp
55,7 → 55,7
Element *theElement = 0;
 
int numArgs = OPS_GetNumRemainingInputArgs();
if (numArgs < 3) {
if (numArgs < 9) {
opserr << "Invalid #args, want: element CompositeElement tag iNode jNode A E I crdTag hinge1 hinge2 \n";
return 0;
}
85,11 → 85,34
UniaxialMaterial *end1 = OPS_getUniaxialMaterial(iData[4]);
UniaxialMaterial *end2 = OPS_getUniaxialMaterial(iData[5]);
 
 
double maxCM = 1e20;
double cmRatio = 0.;
double mass = 0;
while(OPS_GetNumRemainingInputArgs() > 0) {
const char *type = OPS_GetString();
 
if (strcmp(type,"-maxC") == 0) {
if (OPS_GetNumRemainingInputArgs() > 0) {
numData = 1;
if (OPS_GetDoubleInput(&numData,&maxCM) < 0) return 0;
if (OPS_GetNumRemainingInputArgs() > 0) {
numData = 1;
if (OPS_GetDoubleInput(&numData,&cmRatio) < 0) return 0;
}
}
} else if (strcmp(type,"-rho")==0) {
if (OPS_GetNumRemainingInputArgs() > 0) {
numData = 1;
if (OPS_GetDoubleInput(&numData,&mass) < 0) return 0;
}
}
}
 
// Parsing was successful, allocate the material
theElement = new ComponentElement2d(iData[0], dData[0], dData[1], dData[2],
iData[1], iData[2],
*theTrans, end1, end2);
 
*theTrans, end1, end2, mass, maxCM, cmRatio);
if (theElement == 0) {
opserr << "WARNING could not create element of type ComponentElement2d\n";
return 0;
101,8 → 124,10
 
ComponentElement2d::ComponentElement2d()
:Element(0,ELE_TAG_ComponentElement2d),
A(0.0), E(0.0), I(0.0), rho(0.0),
Q(6), q(3), connectedExternalNodes(2), theCoordTransf(0)
A(0.0), E(0.0), I(0.0), rho(0.0),
Q(6), q(3), connectedExternalNodes(2), theCoordTransf(0),
end1Hinge(0), end2Hinge(0),
kTrial(2,2), R(4), uTrial(4), uCommit(4), kb(3,3), kb0(3,3), init(false)
{
// does nothing
q0[0] = 0.0;
121,12 → 146,13
ComponentElement2d::ComponentElement2d(int tag, double a, double e, double i,
int Nd1, int Nd2, CrdTransf &coordTransf,
UniaxialMaterial *end1, UniaxialMaterial *end2,
double r)
double r, double maxC, double cmRat)
:Element(tag,ELE_TAG_ComponentElement2d),
A(a), E(e), I(i), rho(r),
Q(6), q(3), kb(3,3),
Q(6), q(3),
connectedExternalNodes(2), theCoordTransf(0), end1Hinge(0), end2Hinge(0),
kTrial(2,2), R(4), uTrial(4), uCommit(4), init(false)
kTrial(2,2), R(4), uTrial(4), uCommit(4), kb(3,3), kb0(3,3), init(false),
maxCM(maxC), cmRatio(cmRat)
{
connectedExternalNodes(0) = Nd1;
connectedExternalNodes(1) = Nd2;
247,6 → 273,22
EAoverL = A*E/L; // EA/L
EIoverL2 = 2.0*I*E/L; // 2EI/L
EIoverL4 = 4.0*E*I/L; // 4EI/L
 
double k1 = 0.;
if (end1Hinge != 0)
k1 = end1Hinge->getInitialTangent();
double k2 = 0.;
if (end2Hinge != 0)
k2 = end2Hinge->getInitialTangent();
double delta = 1.0/((k1+EIoverL4)*(k2+EIoverL4)-EIoverL2*EIoverL2);
// compute new condensed matrix
kb0(0,0) = EAoverL;
kb0(1,1) = k1 - delta*(k1*k1*(k2+EIoverL4));
kb0(2,2) = k2 - delta*(k2*k2*(k1+EIoverL4));
kb0(1,2) = delta*(k1*k2*EIoverL2);
kb0(2,1) = delta*(k1*k2*EIoverL2);
}
 
int
352,7 → 394,7
bool converged = false;
int count = 0;
int maxCount = 10;
double tol = 1.0e-10;
double tol = 1.0e-14;
 
// iterate, at least once, to remove internal node unbalance
while (converged == false) {
488,9 → 530,9
// determine q = kv + q0
static Vector R(6);
 
q(0) += q0[0];
q(1) += q0[1];
q(2) += q0[2];
// q(0) += q0[0];
// q(1) += q0[1];
// q(2) += q0[2];
kb(0,0) = EAoverL;
kb(1,1) = kTrial(0,0);
504,26 → 546,34
const Matrix &
ComponentElement2d::getInitialStiff(void)
{
double L = theCoordTransf->getInitialLength();
return theCoordTransf->getInitialGlobalStiffMatrix(kb0);
}
 
double k1 = 0.;
if (end1Hinge != 0)
k1 = end1Hinge->getInitialTangent();
double k2 = 0.;
if (end2Hinge != 0)
k2 = end2Hinge->getInitialTangent();
 
double delta = 1.0/((k1+EIoverL4)*(k2+EIoverL4)-EIoverL2*EIoverL2);
const Matrix &
ComponentElement2d::getDamp(void)
{
static Matrix damp(6,6);
damp.Zero();
 
// static Vector theVectorD(6);
// theVectorD = this->getRayleighDampingForces();
// compute new condensed matrix
static Matrix kb0(3,3);
kb0(0,0) = EAoverL;
kb0(1,1) = k1 - delta*(k1*k1*(k2+EIoverL4));
kb0(2,2) = k2 - delta*(k2*k2*(k1+EIoverL4));
kb0(1,2) = delta*(k1*k2*EIoverL2);
kb0(2,1) = delta*(k1*k2*EIoverL2);
return theCoordTransf->getInitialGlobalStiffMatrix(kb0);
// check moment values, if greater than allowed recompute
// if (fabs(theVectorD(2)) > maxCM || fabs(theVectorD(5)) > maxCM)
// return damp;
 
if (alphaM != 0.0)
damp.addMatrix(0.0, this->getMass(), alphaM);
if (betaK != 0.0)
damp.addMatrix(1.0, this->getTangentStiff(), betaK);
if (betaK0 != 0.0)
damp.addMatrix(1.0, this->getInitialStiff(), betaK0);
if (betaKc != 0.0)
damp.addMatrix(1.0, *Kc, betaKc);
 
// return the computed matrix
return damp;
}
 
const Matrix &
661,8 → 711,39
P.addVector(1.0, Q, -1.0);
// add the damping forces if rayleigh damping
if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
P.addVector(1.0, this->getRayleighDampingForces(), 1.0);
if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0) {
static Vector theVectorD(6);
theVectorD = this->getRayleighDampingForces();
// check moment values, if greater than allowed recompute
if (fabs(theVectorD(2)) > maxCM || fabs(theVectorD(5)) > maxCM) {
// opserr << "maxCM: " << maxCM << "\nBEFORE: " << theVectorD;
const Vector &basicVel = theCoordTransf->getBasicTrialVel();
 
kb(0,0) = EAoverL;
kb(1,1) = kTrial(0,0);
kb(2,2) = kTrial(1,1);
kb(1,2) = kTrial(0,1);
kb(2,1) = kTrial(1,0);
Vector localCv(3);
static Vector c0Vec(3);
 
localCv = (betaK + betaKc) * kb * basicVel;
localCv += (betaK0) * kb0 * basicVel;
 
 
if (localCv(1) > maxCM) localCv(1) = maxCM + cmRatio*(localCv(1)-maxCM);
if (localCv(2) > maxCM) localCv(2) = maxCM + cmRatio*(localCv(2)-maxCM);
if (localCv(1) < -maxCM) localCv(1) = -maxCM + cmRatio*(localCv(1)+maxCM);
if (localCv(2) < -maxCM) localCv(2) = -maxCM + cmRatio*(localCv(2)+maxCM);
 
theVectorD = theCoordTransf->getGlobalResistingForce(localCv, c0Vec);
// opserr << this->getTag() << " AFTER: " << theVectorD;
}
P.addVector(1.0, theVectorD, 1.0);
}
// P.addVector(1.0, this->getRayleighDampingForces(), 1.0);
if (rho == 0.0)
return P;
1000,7 → 1081,19
theResponse = new ElementResponse(this, 6, Vector(2));
}
 
else if (strcmp(argv[0],"rayleighForce") == 0 || strcmp(argv[0],"rayleighForces") == 0) {
 
 
output.tag("ResponseType","Px_1");
output.tag("ResponseType","Py_1");
output.tag("ResponseType","Mz_1");
output.tag("ResponseType","Px_2");
output.tag("ResponseType","Py_2");
output.tag("ResponseType","Mz_2");
 
theResponse = new ElementResponse(this, 7, P);
}
 
output.endTag(); // ElementOutput
 
return theResponse;
1063,7 → 1156,12
return eleInfo.setVector(vect2);
 
 
case 7: // basic forces
 
// return eleInfo.setVector(this->getRayleighDampingForces());
return eleInfo.setVector(this->getResistingForceIncInertia());
 
 
default:
return -1;
}
/trunk/SRC/element/componentElement/ComponentElement2d.h
47,7 → 47,7
ComponentElement2d(int tag, double A, double E, double I,
int Nd1, int Nd2, CrdTransf &theTransf,
UniaxialMaterial *end1, UniaxialMaterial *end2,
double rho = 0.0);
double rho = 0.0, double maxCM = 1e20, double cmRatio = 0.);
~ComponentElement2d();
 
const char *getClassType(void) const {return "ComponentElement2d";};
66,6 → 66,7
int update(void);
const Matrix &getTangentStiff(void);
const Matrix &getInitialStiff(void);
const Matrix &getDamp(void);
const Matrix &getMass(void);
 
void zeroLoad(void);
111,6 → 112,7
Vector rTrial;
Vector rCommit;
Matrix kb;
Matrix kb0;
 
static Vector P;
static Matrix K;
120,6 → 122,9
double EAoverL; // EA/L
double EIoverL2; // 2EI/L
double EIoverL4; // 4EI/L
 
double maxCM;
double cmRatio;
};
 
#endif
/trunk/SRC/element/componentElement/Makefile
1,6 → 1,6
include ../../../Makefile.def
 
OBJS = ComponentElement2d.o
OBJS = ComponentElement2d.o ComponentElementDamp2d.o
 
# Compilation control