Skip to content

Commit 3b5c70c

Browse files
mhscottmhscott
authored
Merge pull request OpenSees#1713 from mhscott/mixed-beam-column-matvec
Using addMatrixVector, etc. for mixed beam-column state determination
2 parents 7123c1b + 3720f0d commit 3b5c70c

2 files changed

Lines changed: 224 additions & 72 deletions

File tree

‎SRC/element/mixedBeamColumn/MixedBeamColumn2d.cpp‎

Lines changed: 91 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -862,8 +862,10 @@ int MixedBeamColumn2d::update() {
862862
naturalDisp.Zero();
863863
naturalDisp = crdTransf->getBasicTrialDisp();
864864

865-
naturalIncrDeltaDisp.Zero();
866-
naturalIncrDeltaDisp = naturalDisp - lastNaturalDisp;
865+
//naturalIncrDeltaDisp.Zero();
866+
//naturalIncrDeltaDisp = naturalDisp - lastNaturalDisp;
867+
naturalIncrDeltaDisp = naturalDisp;
868+
naturalIncrDeltaDisp.addVector(1.0, lastNaturalDisp, -1.0);
867869
lastNaturalDisp = naturalDisp;
868870

869871
// Get the numerical integration weights
@@ -883,30 +885,44 @@ int MixedBeamColumn2d::update() {
883885

884886
for( j = 0; j < NDM_SECTION; j++ ){
885887
for( k = 0; k < NDM_NATURAL; k++ ){
886-
nd1T[i](k,j) = nd1[i](j,k);
887-
nd2T[i](k,j) = nd2[i](j,k);
888+
//nd1T[i](k,j) = nd1[i](j,k);
889+
//nd2T[i](k,j) = nd2[i](j,k);
888890
}
889891
}
890892
}
891893

894+
static Vector Vtmp(NDM_NATURAL);
895+
Vtmp = V;
896+
Vtmp.addMatrixVector(1.0, GMH, naturalIncrDeltaDisp, 1.0);
897+
892898
// Update natural force
893899
if (geomLinear) {
894-
naturalForce = naturalForce + Hinv * ( GMH * naturalIncrDeltaDisp + V );
900+
//naturalForce = naturalForce + Hinv * ( GMH * naturalIncrDeltaDisp + V );
901+
naturalForce.addMatrixVector(1.0, Hinv, Vtmp, 1.0);
895902
} else {
896-
naturalForce = naturalForce + Hinv * ( GMH * naturalIncrDeltaDisp + V );
903+
//naturalForce = naturalForce + Hinv * ( GMH * naturalIncrDeltaDisp + V );
904+
naturalForce.addMatrixVector(1.0, Hinv, Vtmp, 1.0);
897905
}
898906

907+
static Vector sectionTmpForce(NDM_SECTION);
908+
899909
// Update sections
900910
for ( i = 0; i < numSections; i++){
901911
// Compute section deformations
902912
sectionForceShapeFcn[i] = nd1[i] * naturalForce;
913+
//sectionForceShapeFcn[i].addMatrixVector(0.0, nd1[i], naturalForce, 1.0);
903914
if (sp != 0) {
904915
const Matrix &s_p = *sp;
905916
for ( j = 0; j < NDM_SECTION; j++ ) {
906917
sectionForceShapeFcn[i](j) += s_p(j,i);
907918
}
908919
}
909-
sectionDefFibers[i] = sectionDefFibers[i] + sectionFlexibility[i] * ( sectionForceShapeFcn[i] - sectionForceFibers[i] );
920+
921+
sectionTmpForce = sectionForceShapeFcn[i];
922+
sectionTmpForce.addVector(1.0, sectionForceFibers[i], -1.0);
923+
924+
//sectionDefFibers[i] = sectionDefFibers[i] + sectionFlexibility[i] * ( sectionForceShapeFcn[i] - sectionForceFibers[i] );
925+
sectionDefFibers[i].addMatrixVector(1.0, sectionFlexibility[i], sectionTmpForce, 1.0);
910926

911927
// Send section deformation to section object
912928
setSectionDeformation(i,sectionDefFibers[i]);
@@ -933,14 +949,46 @@ int MixedBeamColumn2d::update() {
933949
Md.Zero();
934950
Kg.Zero();
935951

952+
static Vector sectionTmpDef(NDM_SECTION);
953+
936954
for( i = 0; i < numSections; i++ ){
937-
V = V + initialLength * wt[i] * nd1T[i] * (sectionDefShapeFcn[i] - sectionDefFibers[i] - sectionFlexibility[i] * ( sectionForceShapeFcn[i] - sectionForceFibers[i] ) );
938-
V2 = V2 + initialLength * wt[i] * nd2T[i] * (sectionDefShapeFcn[i] - sectionDefFibers[i]);
939-
G = G + initialLength * wt[i] * nd1T[i] * nldhat[i];
940-
G2 = G2 + initialLength * wt[i] * nd2T[i] * nldhat[i];
941-
H = H + initialLength * wt[i] * nd1T[i] * sectionFlexibility[i] * nd1[i];
942-
H12 = H12 + initialLength * wt[i] * nd1T[i] * sectionFlexibility[i] * nd2[i];
943-
H22 = H22 + initialLength * wt[i] * nd2T[i] * sectionFlexibility[i] * nd2[i];
955+
sectionTmpDef = sectionDefShapeFcn[i];
956+
sectionTmpDef.addVector(1.0, sectionDefFibers[i], -1.0);
957+
958+
sectionTmpForce = sectionForceShapeFcn[i];
959+
sectionTmpForce.addVector(1.0, sectionForceFibers[i], -1.0);
960+
961+
if (!geomLinear) {
962+
//V2 = V2 + initialLength * wt[i] * nd2T[i] * (sectionDefShapeFcn[i] - sectionDefFibers[i]);
963+
V2.addMatrixTransposeVector(1.0, nd2[i], sectionTmpDef, initialLength*wt[i]);
964+
}
965+
966+
//V = V + initialLength * wt[i] * nd1T[i] * (sectionDefShapeFcn[i] - sectionDefFibers[i] - sectionFlexibility[i] * ( sectionForceShapeFcn[i] - sectionForceFibers[i] ) );
967+
V.addMatrixTransposeVector(1.0, nd1[i], sectionTmpDef, initialLength*wt[i]);
968+
sectionTmpDef.addMatrixVector(0.0, sectionFlexibility[i], sectionTmpForce, 1.0);
969+
V.addMatrixTransposeVector(1.0, nd1[i], sectionTmpDef, -initialLength*wt[i]);
970+
971+
//G = G + initialLength * wt[i] * nd1T[i] * nldhat[i];
972+
G.addMatrixTransposeProduct(1.0, nd1[i], nldhat[i], initialLength*wt[i]);
973+
974+
if (!geomLinear) {
975+
//G2 = G2 + initialLength * wt[i] * nd2T[i] * nldhat[i];
976+
G2.addMatrixTransposeProduct(1.0, nd2[i], nldhat[i], initialLength*wt[i]);
977+
}
978+
979+
//H = H + initialLength * wt[i] * nd1T[i] * sectionFlexibility[i] * nd1[i];
980+
H.addMatrixTripleProduct(1.0, nd1[i], sectionFlexibility[i], initialLength*wt[i]);
981+
982+
if (!geomLinear) {
983+
//H12 = H12 + initialLength * wt[i] * nd1T[i] * sectionFlexibility[i] * nd2[i];
984+
H12.addMatrixTripleProduct(1.0, nd1[i], sectionFlexibility[i], nd2[i], initialLength*wt[i]);
985+
}
986+
987+
if (!geomLinear) {
988+
// H22 = H22 + initialLength * wt[i] * nd2T[i] * sectionFlexibility[i] * nd2[i];
989+
H22.addMatrixTripleProduct(1.0, nd2[i], sectionFlexibility[i], initialLength*wt[i]);
990+
}
991+
944992
if (!geomLinear) {
945993
Kg = Kg + initialLength * wt[i] * this->getKg(i, sectionForceFibers[i](0), currentLength);
946994
// sectionForceFibers[i](0) is the axial load, P
@@ -952,32 +1000,53 @@ int MixedBeamColumn2d::update() {
9521000
invertMatrix(NDM_NATURAL, H, Hinv);
9531001

9541002
// Compute the GMH matrix ( G + Md - H12 ) and its transpose
955-
GMH = G + Md - H12;
1003+
//GMH = G + Md - H12;
1004+
GMH = G;
1005+
if (!geomLinear) {
1006+
GMH.addMatrix(1.0, Md, 1.0);
1007+
GMH.addMatrix(1.0, H12, -1.0);
1008+
}
9561009
//GMH = G; // Omit P-small delta
9571010

9581011
// Compute the transposes of the following matrices: G, G2, GMH
9591012
for( i = 0; i < NDM_NATURAL; i++ ) {
9601013
for( j = 0; j < NDM_NATURAL; j++ ) {
961-
GT(i,j) = G(j,i);
962-
G2T(i,j) = G2(j,i);
963-
GMHT(i,j) = GMH(j,i);
1014+
//GT(i,j) = G(j,i);
1015+
//G2T(i,j) = G2(j,i);
1016+
//GMHT(i,j) = GMH(j,i);
9641017
}
9651018
}
9661019

9671020

9681021
// Define the internal force
9691022
if (geomLinear) {
970-
internalForce = GT * naturalForce + V2 + GMHT * Hinv * V;
1023+
//internalForce = GT * naturalForce + V2 + GMHT * Hinv * V;
1024+
internalForce.addMatrixTransposeVector(0.0, G, naturalForce, 1.0);
1025+
//internalForce.addVector(1.0, V2, 1.0); // Zero if geom linear
1026+
internalForce.addMatrixTransposeVector(1.0, GMH, Hinv * V, 1.0);
9711027
} else {
972-
internalForce = GT * naturalForce + V2 + GMHT * Hinv * V;
1028+
//internalForce = GT * naturalForce + V2 + GMHT * Hinv * V;
1029+
internalForce.addMatrixTransposeVector(0.0, G, naturalForce, 1.0);
1030+
internalForce.addVector(1.0, V2, 1.0);
1031+
internalForce.addMatrixTransposeVector(1.0, GMH, Hinv * V, 1.0);
9731032
}
9741033

9751034

9761035
// Compute the stiffness matrix without the torsion term
9771036
if (geomLinear) {
978-
kv = ( Kg + G2 + G2T - H22 ) + GMHT * Hinv * GMH;
1037+
//kv = ( Kg + G2 + G2T - H22 ) + GMHT * Hinv * GMH;
1038+
kv = Kg;
1039+
//kv.addMatrix(1.0, G2, 1.0); // G2, H22 are zero if geom linear
1040+
//kv.addMatrixTranspose(1.0, G2, 1.0);
1041+
//kv.addMatrixTranspose(1.0, H22, -1.0);
1042+
kv.addMatrixTripleProduct(1.0, GMH, Hinv, 1.0);
9791043
} else {
980-
kv = ( Kg + G2 + G2T - H22 ) + GMHT * Hinv * GMH;
1044+
//kv = ( Kg + G2 + G2T - H22 ) + GMHT * Hinv * GMH;
1045+
kv = Kg;
1046+
kv.addMatrix(1.0, G2, 1.0);
1047+
kv.addMatrixTranspose(1.0, G2, 1.0);
1048+
kv.addMatrixTranspose(1.0, H22, -1.0);
1049+
kv.addMatrixTripleProduct(1.0, GMH, Hinv, 1.0);
9811050
}
9821051

9831052
return 0;

0 commit comments

Comments
 (0)