Rev 4659 | Rev 4968 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 4654 | parduino | 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 | // Created: Chris McGann, UW, 10.2011 |
||
| 22 | // |
||
| 23 | // Description: This file contains the implementation of the SSPbrick class |
||
| 24 | |||
| 25 | #include "SSPbrick.h" |
||
| 26 | |||
| 27 | #include <elementAPI.h> |
||
| 28 | #include <Information.h> |
||
| 29 | #include <ElementResponse.h> |
||
| 30 | #include <ElementalLoad.h> |
||
| 31 | #include <ID.h> |
||
| 32 | #include <Domain.h> |
||
| 33 | #include <Node.h> |
||
| 34 | #include <Channel.h> |
||
| 35 | #include <Message.h> |
||
| 36 | #include <FEM_ObjectBroker.h> |
||
| 37 | #include <Renderer.h> |
||
| 38 | #include <G3Globals.h> |
||
| 39 | #include <ErrorHandler.h> |
||
| 40 | #include <NDMaterial.h> |
||
| 41 | #include <Parameter.h> |
||
| 42 | |||
| 43 | #include <math.h> |
||
| 44 | #include <stdlib.h> |
||
| 45 | #include <stdio.h> |
||
| 46 | |||
| 47 | #define OPS_Export |
||
| 48 | |||
| 49 | static int num_SSPbrick = 0; |
||
| 50 | |||
| 51 | OPS_Export void * |
||
| 52 | OPS_SSPbrick(void) |
||
| 53 | { |
||
| 54 | if (num_SSPbrick == 0) { |
||
| 55 | num_SSPbrick++; |
||
| 56 | OPS_Error("SSPbrick element - Written: C.McGann, P.Arduino, P.Mackenzie-Helnwein, U.Washington\n", 1); |
||
| 57 | } |
||
| 58 | |||
| 59 | // Pointer to an element that will be returned |
||
| 60 | Element *theElement = 0; |
||
| 61 | |||
| 62 | int numRemainingInputArgs = OPS_GetNumRemainingInputArgs(); |
||
| 63 | |||
| 64 | if (numRemainingInputArgs < 10) { |
||
| 65 | opserr << "Invalid #args, want: element SSPbrick eleTag? iNode? jNode? kNode? lNode? mNode? nNode? pNode? qNode? matTag? <b1? b2? b3?>\n"; |
||
| 66 | return 0; |
||
| 67 | } |
||
| 68 | |||
| 69 | int iData[10]; |
||
| 70 | double dData[3]; |
||
| 71 | dData[0] = 0.0; |
||
| 72 | dData[1] = 0.0; |
||
| 73 | dData[2] = 0.0; |
||
| 74 | |||
| 75 | int numData = 10; |
||
| 76 | if (OPS_GetIntInput(&numData, iData) != 0) { |
||
| 77 | opserr << "WARNING invalid integer data: element SSPbrick " << iData[0] << endln; |
||
| 78 | return 0; |
||
| 79 | } |
||
| 80 | |||
| 81 | int matID = iData[9]; |
||
| 82 | NDMaterial *theMaterial = OPS_GetNDMaterial(matID); |
||
| 83 | if (theMaterial == 0) { |
||
| 84 | opserr << "WARNING element SSPbrick " << iData[0] << endln; |
||
| 85 | opserr << " Material: " << matID << "not found\n"; |
||
| 86 | return 0; |
||
| 87 | } |
||
| 88 | |||
| 89 | if (numRemainingInputArgs == 13) { |
||
| 90 | numData = 3; |
||
| 91 | if (OPS_GetDoubleInput(&numData, dData) != 0) { |
||
| 92 | opserr << "WARNING invalid optional data: element SSPbrick " << iData[0] << endln; |
||
| 93 | return 0; |
||
| 94 | } |
||
| 95 | } |
||
| 96 | |||
| 97 | // parsing was successful, allocate the element |
||
| 98 | theElement = new SSPbrick(iData[0], iData[1], iData[2], iData[3], iData[4], iData[5], iData[6], iData[7], iData[8], |
||
| 99 | *theMaterial, dData[0], dData[1], dData[2]); |
||
| 100 | |||
| 101 | if (theElement == 0) { |
||
| 102 | opserr << "WARNING could not create element of type SSPbrick\n"; |
||
| 103 | return 0; |
||
| 104 | } |
||
| 105 | |||
| 106 | return theElement; |
||
| 107 | } |
||
| 108 | |||
| 109 | // full constructor |
||
| 110 | SSPbrick::SSPbrick(int tag, int Nd1, int Nd2, int Nd3, int Nd4, int Nd5, int Nd6, int Nd7, int Nd8, |
||
| 111 | NDMaterial &theMat, double b1, double b2, double b3) |
||
| 112 | :Element(tag,ELE_TAG_SSPbrick), |
||
| 113 | theMaterial(0), |
||
| 114 | mExternalNodes(SSPB_NUM_NODE), |
||
| 115 | mTangentStiffness(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 116 | mInternalForces(SSPB_NUM_DOF), |
||
| 117 | Q(SSPB_NUM_DOF), |
||
| 118 | mMass(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 119 | mNodeCrd(SSPB_NUM_DIM,SSPB_NUM_NODE), |
||
| 120 | mVol(0), |
||
| 121 | Bnot(6,SSPB_NUM_DOF), |
||
| 122 | Kstab(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 123 | xi(8), |
||
| 124 | et(8), |
||
| 125 | ze(8), |
||
| 126 | hut(8), |
||
| 127 | hus(8), |
||
| 128 | hst(8), |
||
| 129 | hstu(8), |
||
| 130 | applyLoad(0) |
||
| 131 | { |
||
| 132 | mExternalNodes(0) = Nd1; |
||
| 133 | mExternalNodes(1) = Nd2; |
||
| 134 | mExternalNodes(2) = Nd3; |
||
| 135 | mExternalNodes(3) = Nd4; |
||
| 136 | mExternalNodes(4) = Nd5; |
||
| 137 | mExternalNodes(5) = Nd6; |
||
| 138 | mExternalNodes(6) = Nd7; |
||
| 139 | mExternalNodes(7) = Nd8; |
||
| 140 | |||
| 141 | b[0] = b1; |
||
| 142 | b[1] = b2; |
||
| 143 | b[2] = b3; |
||
| 144 | |||
| 145 | appliedB[0] = 0.0; |
||
| 146 | appliedB[1] = 0.0; |
||
| 147 | appliedB[2] = 0.0; |
||
| 148 | |||
| 149 | // get copy of the material object |
||
| 150 | NDMaterial *theMatCopy = theMat.getCopy("ThreeDimensional"); |
||
| 151 | if (theMatCopy != 0) { |
||
| 152 | theMaterial = (NDMaterial *)theMatCopy; |
||
| 153 | } else { |
||
| 154 | opserr << "SSPbrick::SSPbrick - failed to get copy of material model\n";; |
||
| 155 | } |
||
| 156 | |||
| 157 | // check material |
||
| 158 | if (theMaterial == 0) { |
||
| 159 | opserr << "SSPbrick::SSPbrick - failed to allocate material model pointer\n"; |
||
| 160 | exit(-1); |
||
| 161 | } |
||
| 162 | } |
||
| 163 | |||
| 164 | // null constructor |
||
| 165 | SSPbrick::SSPbrick() |
||
| 166 | :Element(0,ELE_TAG_SSPbrick), |
||
| 167 | theMaterial(0), |
||
| 168 | mExternalNodes(SSPB_NUM_NODE), |
||
| 169 | mTangentStiffness(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 170 | mInternalForces(SSPB_NUM_DOF), |
||
| 171 | Q(SSPB_NUM_DOF), |
||
| 172 | mMass(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 173 | mNodeCrd(SSPB_NUM_DIM,SSPB_NUM_NODE), |
||
| 174 | mVol(0), |
||
| 175 | Bnot(6,SSPB_NUM_DOF), |
||
| 176 | Kstab(SSPB_NUM_DOF,SSPB_NUM_DOF), |
||
| 177 | xi(8), |
||
| 178 | et(8), |
||
| 179 | ze(8), |
||
| 180 | hut(8), |
||
| 181 | hus(8), |
||
| 182 | hst(8), |
||
| 183 | hstu(8), |
||
| 184 | applyLoad(0) |
||
| 185 | { |
||
| 186 | } |
||
| 187 | |||
| 188 | // destructor |
||
| 189 | SSPbrick::~SSPbrick() |
||
| 190 | { |
||
| 191 | if (theMaterial != 0) { |
||
| 192 | delete theMaterial; |
||
| 193 | } |
||
| 194 | } |
||
| 195 | |||
| 196 | int |
||
| 197 | SSPbrick::getNumExternalNodes(void) const |
||
| 198 | { |
||
| 199 | return SSPB_NUM_NODE; |
||
| 200 | } |
||
| 201 | |||
| 202 | const ID & |
||
| 203 | SSPbrick::getExternalNodes(void) |
||
| 204 | { |
||
| 205 | return mExternalNodes; |
||
| 206 | } |
||
| 207 | |||
| 208 | Node ** |
||
| 209 | SSPbrick::getNodePtrs(void) |
||
| 210 | { |
||
| 211 | return theNodes; |
||
| 212 | } |
||
| 213 | |||
| 214 | int |
||
| 215 | SSPbrick::getNumDOF(void) |
||
| 216 | { |
||
| 217 | return SSPB_NUM_DOF; |
||
| 218 | } |
||
| 219 | |||
| 220 | void |
||
| 221 | SSPbrick::setDomain(Domain *theDomain) |
||
| 222 | { |
||
| 223 | theNodes[0] = theDomain->getNode(mExternalNodes(0)); |
||
| 224 | theNodes[1] = theDomain->getNode(mExternalNodes(1)); |
||
| 225 | theNodes[2] = theDomain->getNode(mExternalNodes(2)); |
||
| 226 | theNodes[3] = theDomain->getNode(mExternalNodes(3)); |
||
| 227 | theNodes[4] = theDomain->getNode(mExternalNodes(4)); |
||
| 228 | theNodes[5] = theDomain->getNode(mExternalNodes(5)); |
||
| 229 | theNodes[6] = theDomain->getNode(mExternalNodes(6)); |
||
| 230 | theNodes[7] = theDomain->getNode(mExternalNodes(7)); |
||
| 231 | |||
| 232 | for (int i = 0; i < 8; i++) { |
||
| 233 | if (theNodes[i] == 0) { |
||
| 234 | return; // don't go any further - otherwise segmentation fault |
||
| 235 | } |
||
| 236 | } |
||
| 237 | Vector mIcrd_1(3); |
||
| 238 | Vector mIcrd_2(3); |
||
| 239 | Vector mIcrd_3(3); |
||
| 240 | Vector mIcrd_4(3); |
||
| 241 | Vector mIcrd_5(3); |
||
| 242 | Vector mIcrd_6(3); |
||
| 243 | Vector mIcrd_7(3); |
||
| 244 | Vector mIcrd_8(3); |
||
| 245 | |||
| 246 | // initialize coordinate vectors |
||
| 247 | mIcrd_1 = theNodes[0]->getCrds(); |
||
| 248 | mIcrd_2 = theNodes[1]->getCrds(); |
||
| 249 | mIcrd_3 = theNodes[2]->getCrds(); |
||
| 250 | mIcrd_4 = theNodes[3]->getCrds(); |
||
| 251 | mIcrd_5 = theNodes[4]->getCrds(); |
||
| 252 | mIcrd_6 = theNodes[5]->getCrds(); |
||
| 253 | mIcrd_7 = theNodes[6]->getCrds(); |
||
| 254 | mIcrd_8 = theNodes[7]->getCrds(); |
||
| 255 | |||
| 256 | // initialize coordinate matrix |
||
| 257 | mNodeCrd(0,0) = mIcrd_1(0); |
||
| 258 | mNodeCrd(1,0) = mIcrd_1(1); |
||
| 259 | mNodeCrd(2,0) = mIcrd_1(2); |
||
| 260 | mNodeCrd(0,1) = mIcrd_2(0); |
||
| 261 | mNodeCrd(1,1) = mIcrd_2(1); |
||
| 262 | mNodeCrd(2,1) = mIcrd_2(2); |
||
| 263 | mNodeCrd(0,2) = mIcrd_3(0); |
||
| 264 | mNodeCrd(1,2) = mIcrd_3(1); |
||
| 265 | mNodeCrd(2,2) = mIcrd_3(2); |
||
| 266 | mNodeCrd(0,3) = mIcrd_4(0); |
||
| 267 | mNodeCrd(1,3) = mIcrd_4(1); |
||
| 268 | mNodeCrd(2,3) = mIcrd_4(2); |
||
| 269 | mNodeCrd(0,4) = mIcrd_5(0); |
||
| 270 | mNodeCrd(1,4) = mIcrd_5(1); |
||
| 271 | mNodeCrd(2,4) = mIcrd_5(2); |
||
| 272 | mNodeCrd(0,5) = mIcrd_6(0); |
||
| 273 | mNodeCrd(1,5) = mIcrd_6(1); |
||
| 274 | mNodeCrd(2,5) = mIcrd_6(2); |
||
| 275 | mNodeCrd(0,6) = mIcrd_7(0); |
||
| 276 | mNodeCrd(1,6) = mIcrd_7(1); |
||
| 277 | mNodeCrd(2,6) = mIcrd_7(2); |
||
| 278 | mNodeCrd(0,7) = mIcrd_8(0); |
||
| 279 | mNodeCrd(1,7) = mIcrd_8(1); |
||
| 280 | mNodeCrd(2,7) = mIcrd_8(2); |
||
| 281 | |||
| 282 | // establish stabilization terms (based on initial state, only need to compute once) |
||
| 283 | GetStab(); |
||
| 284 | |||
| 285 | // call the base-class method |
||
| 286 | this->DomainComponent::setDomain(theDomain); |
||
| 287 | } |
||
| 288 | |||
| 289 | int |
||
| 290 | SSPbrick::commitState(void) |
||
| 291 | { |
||
| 292 | int retVal = 0; |
||
| 293 | // call element commitState to do any base class stuff |
||
| 294 | if ((retVal = this->Element::commitState()) != 0) { |
||
| 295 | opserr << "SSPbrick::commitState() - failed in base class\n"; |
||
| 296 | } |
||
| 297 | retVal = theMaterial->commitState(); |
||
| 298 | |||
| 299 | return retVal; |
||
| 300 | } |
||
| 301 | |||
| 302 | int |
||
| 303 | SSPbrick::revertToLastCommit(void) |
||
| 304 | { |
||
| 305 | return theMaterial->revertToLastCommit(); |
||
| 306 | } |
||
| 307 | |||
| 308 | int |
||
| 309 | SSPbrick::revertToStart(void) |
||
| 310 | { |
||
| 311 | return theMaterial->revertToStart(); |
||
| 312 | } |
||
| 313 | |||
| 314 | int |
||
| 315 | SSPbrick::update(void) |
||
| 316 | // this function updates variables for an incremental step n to n+1 |
||
| 317 | { |
||
| 318 | // get trial displacement |
||
| 319 | const Vector &mDisp_1 = theNodes[0]->getTrialDisp(); |
||
| 320 | const Vector &mDisp_2 = theNodes[1]->getTrialDisp(); |
||
| 321 | const Vector &mDisp_3 = theNodes[2]->getTrialDisp(); |
||
| 322 | const Vector &mDisp_4 = theNodes[3]->getTrialDisp(); |
||
| 323 | const Vector &mDisp_5 = theNodes[4]->getTrialDisp(); |
||
| 324 | const Vector &mDisp_6 = theNodes[5]->getTrialDisp(); |
||
| 325 | const Vector &mDisp_7 = theNodes[6]->getTrialDisp(); |
||
| 326 | const Vector &mDisp_8 = theNodes[7]->getTrialDisp(); |
||
| 327 | |||
| 328 | // assemble displacement vector |
||
| 329 | Vector u(24); |
||
| 330 | u(0) = mDisp_1(0); |
||
| 331 | u(1) = mDisp_1(1); |
||
| 332 | u(2) = mDisp_1(2); |
||
| 333 | u(3) = mDisp_2(0); |
||
| 334 | u(4) = mDisp_2(1); |
||
| 335 | u(5) = mDisp_2(2); |
||
| 336 | u(6) = mDisp_3(0); |
||
| 337 | u(7) = mDisp_3(1); |
||
| 338 | u(8) = mDisp_3(2); |
||
| 339 | u(9) = mDisp_4(0); |
||
| 340 | u(10) = mDisp_4(1); |
||
| 341 | u(11) = mDisp_4(2); |
||
| 342 | u(12) = mDisp_5(0); |
||
| 343 | u(13) = mDisp_5(1); |
||
| 344 | u(14) = mDisp_5(2); |
||
| 345 | u(15) = mDisp_6(0); |
||
| 346 | u(16) = mDisp_6(1); |
||
| 347 | u(17) = mDisp_6(2); |
||
| 348 | u(18) = mDisp_7(0); |
||
| 349 | u(19) = mDisp_7(1); |
||
| 350 | u(20) = mDisp_7(2); |
||
| 351 | u(21) = mDisp_8(0); |
||
| 352 | u(22) = mDisp_8(1); |
||
| 353 | u(23) = mDisp_8(2); |
||
| 354 | |||
| 355 | // compute strain and send it to the material |
||
| 356 | Vector strain(6); |
||
| 357 | strain = Bnot*u; |
||
| 358 | theMaterial->setTrialStrain(strain); |
||
| 359 | |||
| 360 | return 0; |
||
| 361 | } |
||
| 362 | |||
| 363 | const Matrix & |
||
| 364 | SSPbrick::getTangentStiff(void) |
||
| 365 | // this function computes the tangent stiffness matrix for the element |
||
| 366 | { |
||
| 367 | // get material tangent |
||
| 368 | const Matrix &Cmat = theMaterial->getTangent(); |
||
| 369 | |||
| 370 | // full element stiffness matrix |
||
| 371 | mTangentStiffness = Kstab; |
||
| 372 | mTangentStiffness.addMatrixTripleProduct(1.0, Bnot, Cmat, mVol); |
||
| 373 | |||
| 374 | return mTangentStiffness; |
||
| 375 | } |
||
| 376 | |||
| 377 | const Matrix & |
||
| 378 | SSPbrick::getInitialStiff(void) |
||
| 379 | // this function computes the initial tangent stiffness matrix for the element |
||
| 380 | { |
||
| 381 | return getTangentStiff(); |
||
| 382 | } |
||
| 383 | |||
| 384 | const Matrix & |
||
| 385 | SSPbrick::getMass(void) |
||
| 386 | { |
||
| 387 | mMass.Zero(); |
||
| 388 | |||
| 389 | // get mass density from the material |
||
| 390 | double density = theMaterial->getRho(); |
||
| 391 | |||
| 392 | // return zero matrix if density is zero |
||
| 393 | if (density == 0.0) { |
||
| 394 | return mMass; |
||
| 395 | } |
||
| 396 | |||
| 397 | // use jacobian determinant to get nodal mass values |
||
| 398 | double massTerm; |
||
| 399 | for (int i = 0; i < 8; i++) { |
||
| 4925 | parduino | 400 | massTerm = density*J[0]*(1.0 + (J[1]*xi(i) + J[2]*et(i) + J[3]*ze(i) + J[7] + J[8] + J[9])/3.0 |
| 4654 | parduino | 401 | + (J[4]*hut(i) + J[5]*hus(i) + J[6]*hst(i) + J[10]*ze(i) + J[11]*et(i) + J[12]*xi(i) + J[13]*ze(i) + J[14]*et(i) + J[15]*xi(i))/9.0 |
| 402 | + (J[16]*hstu(i) + J[17]*hut(i) + J[18]*hus(i) + J[19]*hst(i))/27.0); |
||
| 403 | mMass(3*i,3*i) += massTerm; |
||
| 404 | mMass(3*i+1,3*i+1) += massTerm; |
||
| 405 | mMass(3*i+2,3*i+2) += massTerm; |
||
| 406 | } |
||
| 407 | |||
| 408 | return mMass; |
||
| 409 | } |
||
| 410 | |||
| 411 | void |
||
| 412 | SSPbrick::zeroLoad(void) |
||
| 413 | { |
||
| 414 | applyLoad = 0; |
||
| 415 | appliedB[0] = 0.0; |
||
| 416 | appliedB[1] = 0.0; |
||
| 417 | appliedB[2] = 0.0; |
||
| 418 | |||
| 419 | return; |
||
| 420 | } |
||
| 421 | |||
| 422 | int |
||
| 423 | SSPbrick::addLoad(ElementalLoad *theLoad, double loadFactor) |
||
| 424 | { |
||
| 425 | // body forces can be applied in a load pattern |
||
| 426 | int type; |
||
| 427 | const Vector &data = theLoad->getData(type, loadFactor); |
||
| 428 | |||
| 429 | if (type == LOAD_TAG_SelfWeight) { |
||
| 430 | applyLoad = 1; |
||
| 431 | appliedB[0] += loadFactor*b[0]; |
||
| 432 | appliedB[1] += loadFactor*b[1]; |
||
| 433 | appliedB[2] += loadFactor*b[2]; |
||
| 434 | return 0; |
||
| 435 | } else { |
||
| 436 | opserr << "SSPbrick::addLoad - load type unknown for ele with tag: " << this->getTag() << endln; |
||
| 437 | return -1; |
||
| 438 | } |
||
| 439 | |||
| 440 | return -1; |
||
| 441 | } |
||
| 442 | |||
| 443 | int |
||
| 444 | SSPbrick::addInertiaLoadToUnbalance(const Vector &accel) |
||
| 445 | { |
||
| 446 | // get mass density from the material |
||
| 447 | double density = theMaterial->getRho(); |
||
| 448 | |||
| 449 | // do nothing if density is zero |
||
| 450 | if (density == 0.0) { |
||
| 451 | return 0; |
||
| 452 | } |
||
| 453 | |||
| 454 | // Get R * accel from the nodes |
||
| 455 | const Vector &Raccel1 = theNodes[0]->getRV(accel); |
||
| 456 | const Vector &Raccel2 = theNodes[1]->getRV(accel); |
||
| 457 | const Vector &Raccel3 = theNodes[2]->getRV(accel); |
||
| 458 | const Vector &Raccel4 = theNodes[3]->getRV(accel); |
||
| 459 | const Vector &Raccel5 = theNodes[4]->getRV(accel); |
||
| 460 | const Vector &Raccel6 = theNodes[5]->getRV(accel); |
||
| 461 | const Vector &Raccel7 = theNodes[6]->getRV(accel); |
||
| 462 | const Vector &Raccel8 = theNodes[7]->getRV(accel); |
||
| 463 | |||
| 464 | static double ra[24]; |
||
| 465 | ra[0] = Raccel1(0); |
||
| 466 | ra[1] = Raccel1(1); |
||
| 467 | ra[2] = Raccel1(2); |
||
| 468 | ra[3] = Raccel2(0); |
||
| 469 | ra[4] = Raccel2(1); |
||
| 470 | ra[5] = Raccel2(2); |
||
| 471 | ra[6] = Raccel3(0); |
||
| 472 | ra[7] = Raccel3(1); |
||
| 473 | ra[8] = Raccel3(2); |
||
| 474 | ra[9] = Raccel4(0); |
||
| 475 | ra[10] = Raccel4(1); |
||
| 476 | ra[11] = Raccel4(2); |
||
| 477 | ra[12] = Raccel5(0); |
||
| 478 | ra[13] = Raccel5(1); |
||
| 479 | ra[14] = Raccel5(2); |
||
| 480 | ra[15] = Raccel6(0); |
||
| 481 | ra[16] = Raccel6(1); |
||
| 482 | ra[17] = Raccel6(2); |
||
| 483 | ra[18] = Raccel7(0); |
||
| 484 | ra[19] = Raccel7(1); |
||
| 485 | ra[20] = Raccel7(2); |
||
| 486 | ra[21] = Raccel8(0); |
||
| 487 | ra[22] = Raccel8(1); |
||
| 488 | ra[23] = Raccel8(2); |
||
| 489 | |||
| 490 | // compute mass matrix |
||
| 491 | this->getMass(); |
||
| 492 | |||
| 493 | for (int i = 0; i < 24; i++) { |
||
| 494 | Q(i) += -mMass(i,i)*ra[i]; |
||
| 495 | } |
||
| 496 | |||
| 497 | return 0; |
||
| 498 | } |
||
| 499 | |||
| 500 | const Vector & |
||
| 501 | SSPbrick::getResistingForce(void) |
||
| 502 | // this function computes the resisting force vector for the element |
||
| 503 | { |
||
| 504 | // get stress from the material |
||
| 505 | Vector mStress = theMaterial->getStress(); |
||
| 506 | |||
| 507 | // get trial displacement |
||
| 508 | const Vector &mDisp_1 = theNodes[0]->getTrialDisp(); |
||
| 509 | const Vector &mDisp_2 = theNodes[1]->getTrialDisp(); |
||
| 510 | const Vector &mDisp_3 = theNodes[2]->getTrialDisp(); |
||
| 511 | const Vector &mDisp_4 = theNodes[3]->getTrialDisp(); |
||
| 512 | const Vector &mDisp_5 = theNodes[4]->getTrialDisp(); |
||
| 513 | const Vector &mDisp_6 = theNodes[5]->getTrialDisp(); |
||
| 514 | const Vector &mDisp_7 = theNodes[6]->getTrialDisp(); |
||
| 515 | const Vector &mDisp_8 = theNodes[7]->getTrialDisp(); |
||
| 516 | |||
| 517 | // assemble displacement vector |
||
| 518 | Vector d(24); |
||
| 519 | d(0) = mDisp_1(0); |
||
| 520 | d(1) = mDisp_1(1); |
||
| 521 | d(2) = mDisp_1(2); |
||
| 522 | d(3) = mDisp_2(0); |
||
| 523 | d(4) = mDisp_2(1); |
||
| 524 | d(5) = mDisp_2(2); |
||
| 525 | d(6) = mDisp_3(0); |
||
| 526 | d(7) = mDisp_3(1); |
||
| 527 | d(8) = mDisp_3(2); |
||
| 528 | d(9) = mDisp_4(0); |
||
| 529 | d(10) = mDisp_4(1); |
||
| 530 | d(11) = mDisp_4(2); |
||
| 531 | d(12) = mDisp_5(0); |
||
| 532 | d(13) = mDisp_5(1); |
||
| 533 | d(14) = mDisp_5(2); |
||
| 534 | d(15) = mDisp_6(0); |
||
| 535 | d(16) = mDisp_6(1); |
||
| 536 | d(17) = mDisp_6(2); |
||
| 537 | d(18) = mDisp_7(0); |
||
| 538 | d(19) = mDisp_7(1); |
||
| 539 | d(20) = mDisp_7(2); |
||
| 540 | d(21) = mDisp_8(0); |
||
| 541 | d(22) = mDisp_8(1); |
||
| 542 | d(23) = mDisp_8(2); |
||
| 543 | |||
| 544 | // add stabilization force to internal force vector |
||
| 545 | mInternalForces = Kstab*d; |
||
| 546 | |||
| 547 | // add internal force from the stress -> fint = Kstab*d + 8*Jo*Bnot'*stress |
||
| 548 | mInternalForces.addMatrixTransposeVector(1.0, Bnot, mStress, mVol); |
||
| 549 | |||
| 550 | // subtract body forces from internal force vector |
||
| 551 | Vector body(3); |
||
| 552 | if (applyLoad == 0) { |
||
| 553 | double polyJac = 0.0; |
||
| 554 | for (int i = 0; i < 8; i++) { |
||
| 4925 | parduino | 555 | polyJac = J[0]*(1.0 + (J[1]*xi(i) + J[2]*et(i) + J[3]*ze(i) + J[7] + J[8] + J[9])/3.0 |
| 4654 | parduino | 556 | + (J[4]*hut(i) + J[5]*hus(i) + J[6]*hst(i) + J[10]*ze(i) + J[11]*et(i) + J[12]*xi(i) + J[13]*ze(i) + J[14]*et(i) + J[15]*xi(i))/9.0 |
| 4925 | parduino | 557 | + (J[16]*hstu(i) + J[17]*hut(i) + J[18]*hus(i) + J[19]*hst(i))/27.0); |
| 4654 | parduino | 558 | mInternalForces(3*i) -= b[0]*polyJac; |
| 559 | mInternalForces(3*i+1) -= b[1]*polyJac; |
||
| 560 | mInternalForces(3*i+2) -= b[2]*polyJac; |
||
| 561 | } |
||
| 562 | } else { |
||
| 563 | double polyJac = 0.0; |
||
| 564 | for (int i = 0; i < 8; i++) { |
||
| 4925 | parduino | 565 | polyJac = J[0]*(1.0 + (J[1]*xi(i) + J[2]*et(i) + J[3]*ze(i) + J[7] + J[8] + J[9])/3.0 |
| 4654 | parduino | 566 | + (J[4]*hut(i) + J[5]*hus(i) + J[6]*hst(i) + J[10]*ze(i) + J[11]*et(i) + J[12]*xi(i) + J[13]*ze(i) + J[14]*et(i) + J[15]*xi(i))/9.0 |
| 4925 | parduino | 567 | + (J[16]*hstu(i) + J[17]*hut(i) + J[18]*hus(i) + J[19]*hst(i))/27.0); |
| 4654 | parduino | 568 | mInternalForces(3*i) -= appliedB[0]*polyJac; |
| 569 | mInternalForces(3*i+1) -= appliedB[1]*polyJac; |
||
| 570 | mInternalForces(3*i+2) -= appliedB[2]*polyJac; |
||
| 571 | } |
||
| 572 | } |
||
| 573 | |||
| 574 | // inertial unbalance load |
||
| 575 | mInternalForces.addVector(1.0, Q, -1.0); |
||
| 576 | |||
| 577 | return mInternalForces; |
||
| 578 | } |
||
| 579 | |||
| 580 | const Vector & |
||
| 581 | SSPbrick::getResistingForceIncInertia() |
||
| 582 | { |
||
| 583 | // get mass density from the material |
||
| 584 | double density = theMaterial->getRho(); |
||
| 585 | |||
| 586 | // if density is zero only add damping terms |
||
| 587 | if (density == 0.0) { |
||
| 588 | this->getResistingForce(); |
||
| 589 | |||
| 590 | // add the damping forces if rayleigh damping |
||
| 591 | if (betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0) { |
||
| 592 | mInternalForces += this->getRayleighDampingForces(); |
||
| 593 | } |
||
| 594 | |||
| 595 | return mInternalForces; |
||
| 596 | } |
||
| 597 | |||
| 598 | const Vector &accel1 = theNodes[0]->getTrialAccel(); |
||
| 599 | const Vector &accel2 = theNodes[1]->getTrialAccel(); |
||
| 600 | const Vector &accel3 = theNodes[2]->getTrialAccel(); |
||
| 601 | const Vector &accel4 = theNodes[3]->getTrialAccel(); |
||
| 602 | const Vector &accel5 = theNodes[4]->getTrialAccel(); |
||
| 603 | const Vector &accel6 = theNodes[5]->getTrialAccel(); |
||
| 604 | const Vector &accel7 = theNodes[6]->getTrialAccel(); |
||
| 605 | const Vector &accel8 = theNodes[7]->getTrialAccel(); |
||
| 606 | |||
| 607 | static double a[24]; |
||
| 608 | a[0] = accel1(0); |
||
| 609 | a[1] = accel1(1); |
||
| 610 | a[2] = accel1(2); |
||
| 611 | a[3] = accel2(0); |
||
| 612 | a[4] = accel2(1); |
||
| 613 | a[5] = accel2(2); |
||
| 614 | a[6] = accel3(0); |
||
| 615 | a[7] = accel3(1); |
||
| 616 | a[8] = accel3(2); |
||
| 617 | a[9] = accel4(0); |
||
| 618 | a[10] = accel4(1); |
||
| 619 | a[11] = accel4(2); |
||
| 620 | a[12] = accel5(0); |
||
| 621 | a[13] = accel5(1); |
||
| 622 | a[14] = accel5(2); |
||
| 623 | a[15] = accel6(0); |
||
| 624 | a[16] = accel6(1); |
||
| 625 | a[17] = accel6(2); |
||
| 626 | a[18] = accel7(0); |
||
| 627 | a[19] = accel7(1); |
||
| 628 | a[20] = accel7(2); |
||
| 629 | a[21] = accel8(0); |
||
| 630 | a[22] = accel8(1); |
||
| 631 | a[23] = accel8(2); |
||
| 632 | |||
| 633 | // compute current resisting force |
||
| 634 | this->getResistingForce(); |
||
| 635 | |||
| 636 | // compute mass matrix |
||
| 637 | this->getMass(); |
||
| 638 | |||
| 639 | for (int i = 0; i < 8; i++) { |
||
| 640 | mInternalForces(i) += mMass(i,i)*a[i]; |
||
| 641 | } |
||
| 642 | |||
| 643 | // add the damping forces if rayleigh damping |
||
| 644 | if (betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0) { |
||
| 645 | mInternalForces += this->getRayleighDampingForces(); |
||
| 646 | } |
||
| 647 | |||
| 648 | return mInternalForces; |
||
| 649 | } |
||
| 650 | |||
| 651 | int |
||
| 652 | SSPbrick::sendSelf(int commitTag, Channel &theChannel) |
||
| 653 | { |
||
| 654 | int res = 0; |
||
| 655 | |||
| 656 | // note: we don't check for dataTag == 0 for Element |
||
| 657 | // objects as that is taken care of in a commit by the Domain |
||
| 658 | // object - don't want to have to do the check if sending data |
||
| 659 | int dataTag = this->getDbTag(); |
||
| 660 | |||
| 661 | // SSPbrick packs its data into a Vector and sends this to theChannel |
||
| 662 | // along with its dbTag and the commitTag passed in the arguments |
||
| 663 | static Vector data(6); |
||
| 664 | data(0) = this->getTag(); |
||
| 665 | data(1) = b[0]; |
||
| 666 | data(2) = b[1]; |
||
| 667 | data(3) = b[2]; |
||
| 668 | data(4) = theMaterial->getClassTag(); |
||
| 669 | |||
| 670 | // Now quad sends the ids of its materials |
||
| 671 | int matDbTag = theMaterial->getDbTag(); |
||
| 672 | |||
| 673 | static ID idData(12); |
||
| 674 | |||
| 675 | // NOTE: we do have to ensure that the material has a database |
||
| 676 | // tag if we are sending to a database channel. |
||
| 677 | if (matDbTag == 0) { |
||
| 678 | matDbTag = theChannel.getDbTag(); |
||
| 679 | if (matDbTag != 0) |
||
| 680 | theMaterial->setDbTag(matDbTag); |
||
| 681 | } |
||
| 682 | data(5) = matDbTag; |
||
| 683 | |||
| 684 | res += theChannel.sendVector(dataTag, commitTag, data); |
||
| 685 | if (res < 0) { |
||
| 686 | opserr << "WARNING SSPbrick::sendSelf() - " << this->getTag() << " failed to send Vector\n"; |
||
| 687 | return res; |
||
| 688 | } |
||
| 689 | |||
| 690 | // SSPbrick then sends the tags of its four nodes |
||
| 691 | res += theChannel.sendID(dataTag, commitTag, mExternalNodes); |
||
| 692 | if (res < 0) { |
||
| 693 | opserr << "WARNING SSPbrick::sendSelf() - " << this->getTag() << " failed to send ID\n"; |
||
| 694 | return res; |
||
| 695 | } |
||
| 696 | |||
| 697 | // finally, SSPbrick asks its material object to send itself |
||
| 698 | res = theMaterial->sendSelf(commitTag, theChannel); |
||
| 699 | if (res < 0) { |
||
| 700 | opserr << "WARNING SSPbrick::sendSelf() - " << this->getTag() << " failed to send its Material\n"; |
||
| 701 | return -3; |
||
| 702 | } |
||
| 703 | |||
| 704 | return 0; |
||
| 705 | } |
||
| 706 | |||
| 707 | int |
||
| 708 | SSPbrick::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) |
||
| 709 | { |
||
| 710 | int res = 0; |
||
| 711 | int dataTag = this->getDbTag(); |
||
| 712 | |||
| 713 | // SSPbrick creates a Vector, receives the Vector and then sets the |
||
| 714 | // internal data with the data in the Vector |
||
| 715 | static Vector data(6); |
||
| 716 | res += theChannel.recvVector(dataTag, commitTag, data); |
||
| 717 | if (res < 0) { |
||
| 718 | opserr << "WARNING SSPbrick::recvSelf() - failed to receive Vector\n"; |
||
| 719 | return res; |
||
| 720 | } |
||
| 721 | |||
| 722 | this->setTag((int)data(0)); |
||
| 723 | b[0] = data(1); |
||
| 724 | b[1] = data(2); |
||
| 725 | b[2] = data(3); |
||
| 726 | |||
| 727 | |||
| 728 | // SSPbrick now receives the tags of its four external nodes |
||
| 729 | res += theChannel.recvID(dataTag, commitTag, mExternalNodes); |
||
| 730 | if (res < 0) { |
||
| 731 | opserr << "WARNING SSPbrick::recvSelf() - " << this->getTag() << " failed to receive ID\n"; |
||
| 732 | return res; |
||
| 733 | } |
||
| 734 | |||
| 735 | // finally, SSPbrick creates a material object of the correct type, sets its |
||
| 736 | // database tag, and asks this new object to receive itself |
||
| 737 | int matClass = (int)data(4); |
||
| 738 | int matDb = (int)data(5); |
||
| 739 | |||
| 740 | // check if material object exists and that it is the right type |
||
| 741 | if ((theMaterial == 0) || (theMaterial->getClassTag() != matClass)) { |
||
| 742 | |||
| 743 | // if old one, delete it |
||
| 744 | if (theMaterial != 0) |
||
| 745 | delete theMaterial; |
||
| 746 | |||
| 747 | // create new material object |
||
| 748 | NDMaterial *theMatCopy = theBroker.getNewNDMaterial(matClass); |
||
| 749 | theMaterial = (NDMaterial *)theMatCopy; |
||
| 750 | |||
| 751 | if (theMaterial == 0) { |
||
| 752 | opserr << "WARNING SSPbrick::recvSelf() - " << this->getTag() |
||
| 753 | << " failed to get a blank Material of type " << matClass << endln; |
||
| 754 | return -3; |
||
| 755 | } |
||
| 756 | } |
||
| 757 | |||
| 758 | // NOTE: we set the dbTag before we receive the material |
||
| 759 | theMaterial->setDbTag(matDb); |
||
| 760 | res = theMaterial->recvSelf(commitTag, theChannel, theBroker); |
||
| 761 | if (res < 0) { |
||
| 762 | opserr << "WARNING SSPbrick::recvSelf() - " << this->getTag() << " failed to receive its Material\n"; |
||
| 763 | return -3; |
||
| 764 | } |
||
| 765 | |||
| 766 | return 0; |
||
| 767 | } |
||
| 768 | |||
| 769 | int |
||
| 770 | SSPbrick::displaySelf(Renderer &theViewer, int displayMode, float fact) |
||
| 771 | { |
||
| 772 | return 0; |
||
| 773 | } |
||
| 774 | |||
| 775 | void |
||
| 776 | SSPbrick::Print(OPS_Stream &s, int flag) |
||
| 777 | { |
||
| 778 | opserr << "SSPbrick, element id: " << this->getTag() << endln; |
||
| 779 | opserr << " Connected external nodes: "; |
||
| 780 | for (int i = 0; i < SSPB_NUM_NODE; i++) { |
||
| 781 | opserr << mExternalNodes(i) << " "; |
||
| 782 | } |
||
| 783 | return; |
||
| 784 | } |
||
| 785 | |||
| 786 | Response* |
||
| 787 | SSPbrick::setResponse(const char **argv, int argc, OPS_Stream &eleInfo) |
||
| 788 | { |
||
| 789 | // no special recorders for this element, call the method in the material class |
||
| 790 | return theMaterial->setResponse(argv, argc, eleInfo); |
||
| 791 | } |
||
| 792 | |||
| 793 | int |
||
| 794 | SSPbrick::getResponse(int responseID, Information &eleInfo) |
||
| 795 | { |
||
| 796 | // no special recorders for this element, call the method in the material class |
||
| 797 | return theMaterial->getResponse(responseID, eleInfo); |
||
| 798 | } |
||
| 799 | |||
| 800 | int |
||
| 801 | SSPbrick::setParameter(const char **argv, int argc, Parameter ¶m) |
||
| 802 | { |
||
| 803 | if (argc < 1) { |
||
| 804 | return -1; |
||
| 805 | } |
||
| 806 | |||
| 807 | int res = -1; |
||
| 808 | |||
| 809 | if (strcmp(argv[0],"materialState") == 0) { |
||
| 810 | return param.addObject(5,this); |
||
| 811 | } |
||
| 812 | |||
| 813 | // quad pressure loading |
||
| 814 | if (strcmp(argv[0],"pressure") == 0) { |
||
| 815 | return param.addObject(2, this); |
||
| 816 | } |
||
| 817 | |||
| 818 | // a material parameter |
||
| 819 | else if (strstr(argv[0],"material") != 0) { |
||
| 820 | |||
| 821 | if (argc < 3) { |
||
| 822 | return -1; |
||
| 823 | } |
||
| 824 | |||
| 825 | int pointNum = atoi(argv[1]); |
||
| 826 | if (pointNum > 0 && pointNum <= 4) { |
||
| 827 | return theMaterial->setParameter(&argv[2], argc-2, param); |
||
| 828 | } else { |
||
| 829 | return -1; |
||
| 830 | } |
||
| 831 | } |
||
| 832 | |||
| 833 | // otherwise it could be just a forall material parameter |
||
| 834 | else { |
||
| 835 | |||
| 836 | int matRes = res; |
||
| 837 | matRes = theMaterial->setParameter(argv, argc, param); |
||
| 838 | |||
| 839 | if (matRes != -1) { |
||
| 840 | res = matRes; |
||
| 841 | } |
||
| 842 | } |
||
| 843 | |||
| 844 | return res; |
||
| 845 | } |
||
| 846 | |||
| 847 | int |
||
| 848 | SSPbrick::updateParameter(int parameterID, Information &info) |
||
| 849 | { |
||
| 850 | int res = -1; |
||
| 851 | int matRes = res; |
||
| 852 | switch (parameterID) { |
||
| 853 | case -1: |
||
| 854 | return -1; |
||
| 855 | |||
| 856 | case 1: |
||
| 857 | matRes = theMaterial->updateParameter(parameterID, info); |
||
| 858 | if (matRes != -1) { |
||
| 859 | res = matRes; |
||
| 860 | } |
||
| 861 | return res; |
||
| 862 | |||
| 863 | case 2: |
||
| 864 | //pressure = info.theDouble; |
||
| 865 | //this->setPressureLoadAtNodes(); // update consistent nodal loads |
||
| 866 | return 0; |
||
| 867 | |||
| 868 | case 5: |
||
| 869 | matRes = theMaterial->updateParameter(parameterID, info); |
||
| 870 | if (matRes != -1) { |
||
| 871 | res = matRes; |
||
| 872 | } |
||
| 873 | return res; |
||
| 874 | |||
| 875 | default: |
||
| 876 | return -1; |
||
| 877 | } |
||
| 878 | } |
||
| 879 | |||
| 880 | void |
||
| 881 | SSPbrick::GetStab(void) |
||
| 882 | // this function computes the stabilization stiffness matrix for the element |
||
| 883 | { |
||
| 884 | Matrix Mben(12,24); |
||
| 885 | Matrix FCF(12,12); |
||
| 886 | Matrix dNloc(8,3); |
||
| 887 | Matrix dNmod(8,3); |
||
| 888 | Matrix Jmat(3,3); |
||
| 889 | Matrix Jinv(3,3); |
||
| 890 | Matrix G(8,8); |
||
| 891 | Matrix I8(8,8); |
||
| 892 | |||
| 893 | // define local coord and hourglass vectors |
||
| 894 | Vector gst(8); |
||
| 895 | Vector gut(8); |
||
| 896 | Vector gus(8); |
||
| 897 | Vector gstu(8); |
||
| 898 | |||
| 899 | xi(0) = -0.125; |
||
| 900 | xi(1) = 0.125; |
||
| 901 | xi(2) = 0.125; |
||
| 902 | xi(3) = -0.125; |
||
| 903 | xi(4) = -0.125; |
||
| 904 | xi(5) = 0.125; |
||
| 905 | xi(6) = 0.125; |
||
| 906 | xi(7) = -0.125; |
||
| 907 | |||
| 908 | et(0) = -0.125; |
||
| 909 | et(1) = -0.125; |
||
| 910 | et(2) = 0.125; |
||
| 911 | et(3) = 0.125; |
||
| 912 | et(4) = -0.125; |
||
| 913 | et(5) = -0.125; |
||
| 914 | et(6) = 0.125; |
||
| 915 | et(7) = 0.125; |
||
| 916 | |||
| 917 | ze(0) = -0.125; |
||
| 918 | ze(1) = -0.125; |
||
| 919 | ze(2) = -0.125; |
||
| 920 | ze(3) = -0.125; |
||
| 921 | ze(4) = 0.125; |
||
| 922 | ze(5) = 0.125; |
||
| 923 | ze(6) = 0.125; |
||
| 924 | ze(7) = 0.125; |
||
| 925 | |||
| 926 | hst(0) = 0.125; |
||
| 927 | hst(1) = -0.125; |
||
| 928 | hst(2) = 0.125; |
||
| 929 | hst(3) = -0.125; |
||
| 930 | hst(4) = 0.125; |
||
| 931 | hst(5) = -0.125; |
||
| 932 | hst(6) = 0.125; |
||
| 933 | hst(7) = -0.125; |
||
| 934 | |||
| 935 | hut(0) = 0.125; |
||
| 936 | hut(1) = 0.125; |
||
| 937 | hut(2) = -0.125; |
||
| 938 | hut(3) = -0.125; |
||
| 939 | hut(4) = -0.125; |
||
| 940 | hut(5) = -0.125; |
||
| 941 | hut(6) = 0.125; |
||
| 942 | hut(7) = 0.125; |
||
| 943 | |||
| 944 | hus(0) = 0.125; |
||
| 945 | hus(1) = -0.125; |
||
| 946 | hus(2) = -0.125; |
||
| 947 | hus(3) = 0.125; |
||
| 948 | hus(4) = -0.125; |
||
| 949 | hus(5) = 0.125; |
||
| 950 | hus(6) = 0.125; |
||
| 951 | hus(7) = -0.125; |
||
| 952 | |||
| 953 | hstu(0) = -0.125; |
||
| 954 | hstu(1) = 0.125; |
||
| 955 | hstu(2) = -0.125; |
||
| 956 | hstu(3) = 0.125; |
||
| 957 | hstu(4) = 0.125; |
||
| 958 | hstu(5) = -0.125; |
||
| 959 | hstu(6) = 0.125; |
||
| 960 | hstu(7) = -0.125; |
||
| 961 | |||
| 962 | // shape function derivatives (local crd) at center |
||
| 963 | dNloc(0,0) = -0.125; |
||
| 964 | dNloc(1,0) = 0.125; |
||
| 965 | dNloc(2,0) = 0.125; |
||
| 966 | dNloc(3,0) = -0.125; |
||
| 967 | dNloc(4,0) = -0.125; |
||
| 968 | dNloc(5,0) = 0.125; |
||
| 969 | dNloc(6,0) = 0.125; |
||
| 970 | dNloc(7,0) = -0.125; |
||
| 971 | dNloc(0,1) = -0.125; |
||
| 972 | dNloc(1,1) = -0.125; |
||
| 973 | dNloc(2,1) = 0.125; |
||
| 974 | dNloc(3,1) = 0.125; |
||
| 975 | dNloc(4,1) = -0.125; |
||
| 976 | dNloc(5,1) = -0.125; |
||
| 977 | dNloc(6,1) = 0.125; |
||
| 978 | dNloc(7,1) = 0.125; |
||
| 979 | dNloc(0,2) = -0.125; |
||
| 980 | dNloc(1,2) = -0.125; |
||
| 981 | dNloc(2,2) = -0.125; |
||
| 982 | dNloc(3,2) = -0.125; |
||
| 983 | dNloc(4,2) = 0.125; |
||
| 984 | dNloc(5,2) = 0.125; |
||
| 985 | dNloc(6,2) = 0.125; |
||
| 986 | dNloc(7,2) = 0.125; |
||
| 987 | |||
| 988 | // jacobian matrix |
||
| 989 | Jmat = mNodeCrd*dNloc; |
||
| 990 | // inverse of jacobian matrix |
||
| 991 | Jmat.Invert(Jinv); |
||
| 992 | |||
| 993 | // nodal coordinate vectors |
||
| 994 | Vector x(8); |
||
| 995 | Vector y(8); |
||
| 996 | Vector z(8); |
||
| 997 | |||
| 998 | x(0) = mNodeCrd(0,0); |
||
| 999 | x(1) = mNodeCrd(0,1); |
||
| 1000 | x(2) = mNodeCrd(0,2); |
||
| 1001 | x(3) = mNodeCrd(0,3); |
||
| 1002 | x(4) = mNodeCrd(0,4); |
||
| 1003 | x(5) = mNodeCrd(0,5); |
||
| 1004 | x(6) = mNodeCrd(0,6); |
||
| 1005 | x(7) = mNodeCrd(0,7); |
||
| 1006 | |||
| 1007 | y(0) = mNodeCrd(1,0); |
||
| 1008 | y(1) = mNodeCrd(1,1); |
||
| 1009 | y(2) = mNodeCrd(1,2); |
||
| 1010 | y(3) = mNodeCrd(1,3); |
||
| 1011 | y(4) = mNodeCrd(1,4); |
||
| 1012 | y(5) = mNodeCrd(1,5); |
||
| 1013 | y(6) = mNodeCrd(1,6); |
||
| 1014 | y(7) = mNodeCrd(1,7); |
||
| 1015 | |||
| 1016 | z(0) = mNodeCrd(2,0); |
||
| 1017 | z(1) = mNodeCrd(2,1); |
||
| 1018 | z(2) = mNodeCrd(2,2); |
||
| 1019 | z(3) = mNodeCrd(2,3); |
||
| 1020 | z(4) = mNodeCrd(2,4); |
||
| 1021 | z(5) = mNodeCrd(2,5); |
||
| 1022 | z(6) = mNodeCrd(2,6); |
||
| 1023 | z(7) = mNodeCrd(2,7); |
||
| 1024 | |||
| 1025 | // define coefficent terms for jacobian determinant |
||
| 1026 | double a1 = x^xi; |
||
| 1027 | double a2 = x^et; |
||
| 1028 | double a3 = x^ze; |
||
| 1029 | double a4 = x^hut; |
||
| 1030 | double a5 = x^hus; |
||
| 1031 | double a6 = x^hst; |
||
| 1032 | double a7 = x^hstu; |
||
| 1033 | |||
| 1034 | double b1 = y^xi; |
||
| 1035 | double b2 = y^et; |
||
| 1036 | double b3 = y^ze; |
||
| 1037 | double b4 = y^hut; |
||
| 1038 | double b5 = y^hus; |
||
| 1039 | double b6 = y^hst; |
||
| 1040 | double b7 = y^hstu; |
||
| 1041 | |||
| 1042 | double c1 = z^xi; |
||
| 1043 | double c2 = z^et; |
||
| 1044 | double c3 = z^ze; |
||
| 1045 | double c4 = z^hut; |
||
| 1046 | double c5 = z^hus; |
||
| 1047 | double c6 = z^hst; |
||
| 1048 | double c7 = z^hstu; |
||
| 1049 | |||
| 1050 | // define coefficient vectors for jacobian determinant |
||
| 1051 | Vector e1(3); |
||
| 1052 | Vector e2(3); |
||
| 1053 | Vector e3(3); |
||
| 1054 | Vector e4(3); |
||
| 1055 | Vector e5(3); |
||
| 1056 | Vector e6(3); |
||
| 1057 | Vector e7(3); |
||
| 1058 | |||
| 1059 | e1(0) = a1; e1(1) = b1; e1(2) = c1; |
||
| 1060 | e2(0) = a2; e2(1) = b2; e2(2) = c2; |
||
| 1061 | e3(0) = a3; e3(1) = b3; e3(2) = c3; |
||
| 1062 | e4(0) = a4; e4(1) = b4; e4(2) = c4; |
||
| 1063 | e5(0) = a5; e5(1) = b5; e5(2) = c5; |
||
| 1064 | e6(0) = a6; e6(1) = b6; e6(2) = c6; |
||
| 1065 | e7(0) = a7; e7(1) = b7; e7(2) = c7; |
||
| 1066 | |||
| 1067 | // jacobian determinant terms |
||
| 1068 | J[0] = e1^(CrossProduct(e2,e3)); |
||
| 1069 | J[1] = (e1^(CrossProduct(e2,e5))) + (e1^(CrossProduct(e6,e3))); |
||
| 1070 | J[2] = (e1^(CrossProduct(e2,e4))) + (e6^(CrossProduct(e2,e3))); |
||
| 1071 | J[3] = (e5^(CrossProduct(e2,e3))) + (e1^(CrossProduct(e4,e3))); |
||
| 4925 | parduino | 1072 | J[1] = 0.0; |
| 1073 | J[2] = 0.0; |
||
| 1074 | J[3] = 0.0; |
||
| 4654 | parduino | 1075 | J[4] = (e7^(CrossProduct(e2,e3))) + (e4^(CrossProduct(e5,e2))) + (e4^(CrossProduct(e3,e6))); |
| 1076 | J[5] = (e1^(CrossProduct(e7,e3))) + (e4^(CrossProduct(e5,e1))) + (e3^(CrossProduct(e5,e6))); |
||
| 1077 | J[6] = (e1^(CrossProduct(e2,e7))) + (e4^(CrossProduct(e1,e6))) + (e2^(CrossProduct(e5,e6))); |
||
| 1078 | J[7] = -1.0*e1^(CrossProduct(e5,e6)); |
||
| 1079 | J[8] = -1.0*e4^(CrossProduct(e2,e6)); |
||
| 1080 | J[9] = -1.0*e4^(CrossProduct(e5,e3)); |
||
| 1081 | J[10] = e2^(CrossProduct(e4,e7)); |
||
| 1082 | J[11] = -1.0*e3^(CrossProduct(e4,e7)); |
||
| 1083 | J[12] = e3^(CrossProduct(e5,e7)); |
||
| 1084 | J[13] = -1.0*e1^(CrossProduct(e5,e7)); |
||
| 1085 | J[14] = e1^(CrossProduct(e6,e7)); |
||
| 1086 | J[15] = -1.0*e2^(CrossProduct(e6,e7)); |
||
| 1087 | J[16] = 2.0*e4^(CrossProduct(e5,e6)); |
||
| 1088 | J[17] = e7^(CrossProduct(e5,e6)); |
||
| 1089 | J[18] = e4^(CrossProduct(e7,e6)); |
||
| 1090 | J[19] = e4^(CrossProduct(e5,e7)); |
||
| 1091 | |||
| 1092 | // combined jacobian terms |
||
| 1093 | double J0789 = 8.0*(J[0]/3.0 + J[7]/5.0 + J[8]/9.0 + J[9]/9.0); |
||
| 1094 | double J0879 = 8.0*(J[0]/3.0 + J[8]/5.0 + J[7]/9.0 + J[9]/9.0); |
||
| 1095 | double J0978 = 8.0*(J[0]/3.0 + J[9]/5.0 + J[7]/9.0 + J[8]/9.0); |
||
| 1096 | double J417 = 8.0*(J[4]/9.0 + J[17]/27.0); |
||
| 1097 | double J518 = 8.0*(J[5]/9.0 + J[18]/27.0); |
||
| 1098 | double J619 = 8.0*(J[6]/9.0 + J[19]/27.0); |
||
| 1099 | double J11215 = 8.0*(J[1]/9.0 + J[12]/15.0 + J[15]/27.0); |
||
| 1100 | double J11512 = 8.0*(J[1]/9.0 + J[15]/15.0 + J[12]/27.0); |
||
| 1101 | double J21114 = 8.0*(J[2]/9.0 + J[11]/15.0 + J[14]/27.0); |
||
| 1102 | double J21411 = 8.0*(J[2]/9.0 + J[14]/15.0 + J[11]/27.0); |
||
| 1103 | double J31013 = 8.0*(J[3]/9.0 + J[10]/15.0 + J[13]/27.0); |
||
| 1104 | double J31310 = 8.0*(J[3]/9.0 + J[13]/15.0 + J[10]/27.0); |
||
| 4925 | parduino | 1105 | double J789 = 8.0*(J[0]/9.0 + J[7]/15.0 + J[8]/15.0 + J[9]/27.0); |
| 4654 | parduino | 1106 | double J897 = 8.0*(J[0]/9.0 + J[8]/15.0 + J[9]/15.0 + J[7]/27.0); |
| 1107 | double J798 = 8.0*(J[0]/9.0 + J[7]/15.0 + J[9]/15.0 + J[8]/27.0); |
||
| 4925 | parduino | 1108 | double J174 = 8.0*(J[4]/27.0 + 64.0*J[17]/45.0); |
| 1109 | double J185 = 8.0*(J[5]/27.0 + 64.0*J[18]/45.0); |
||
| 1110 | double J196 = 8.0*(J[6]/27.0 + 64.0*J[19]/45.0); |
||
| 1111 | double J16 = 8.0*J[16]/27.0; |
||
| 4654 | parduino | 1112 | |
| 1113 | // compute element volume |
||
| 4925 | parduino | 1114 | mVol = 8.0*(J[0] + (J[7] + J[8] + J[9])/3.0); |
| 4654 | parduino | 1115 | |
| 1116 | // kinematic vectors |
||
| 1117 | Vector bx(8); |
||
| 1118 | Vector by(8); |
||
| 1119 | Vector bz(8); |
||
| 1120 | |||
| 1121 | bx = (8.0*((b2*c3-c2*b3)*xi + (b3*c1-c3*b1)*et + (b1*c2-c1*b2)*ze) + (8.0/3.0)*((b6*c5-c6*b5)*xi + (b4*c6-c4*b6)*et + (b5*c4-c5*b4)*ze |
||
| 1122 | + (b5*c1-c5*b1 + b2*c4-c2*b4)*hst + (b6*c2-c6*b2 + b3*c5-c3*b5)*hut + (b1*c6-c1*b6 + b4*c3-c4*b3)*hus))/mVol; |
||
| 1123 | by = (8.0*((c2*a3-a2*c3)*xi + (c3*a1-a3*c1)*et + (c1*a2-a1*c2)*ze) + (8.0/3.0)*((c6*a5-a6*c5)*xi + (c4*a6-a4*c6)*et + (c5*a4-a5*c4)*ze |
||
| 1124 | + (c5*a1-a5*c1 + c2*a4-a2*c4)*hst + (c6*a2-a6*c2 + c3*a5-a3*c5)*hut + (c1*a6-a1*c6 + c4*a3-a4*c3)*hus))/mVol; |
||
| 1125 | bz = (8.0*((a2*b3-b2*a3)*xi + (a3*b1-b3*a1)*et + (a1*b2-b1*a2)*ze) + (8.0/3.0)*((a6*b5-b6*a5)*xi + (a4*b6-b4*a6)*et + (a5*b4-b5*a4)*ze |
||
| 1126 | + (a5*b1-b5*a1 + a2*b4-b2*a4)*hst + (a6*b2-b6*a2 + a3*b5-b3*a5)*hut + (a1*b6-b1*a6 + a4*b3-b4*a3)*hus))/mVol; |
||
| 1127 | |||
| 1128 | for (int i = 0; i < 8; i++) { |
||
| 1129 | dNmod(i,0) = bx(i); |
||
| 1130 | dNmod(i,1) = by(i); |
||
| 1131 | dNmod(i,2) = bz(i); |
||
| 1132 | } |
||
| 1133 | |||
| 1134 | // compute hourglass transformation matrix G |
||
| 1135 | I8.Zero(); |
||
| 1136 | I8(0,0) = 1.0; |
||
| 1137 | I8(1,1) = 1.0; |
||
| 1138 | I8(2,2) = 1.0; |
||
| 1139 | I8(3,3) = 1.0; |
||
| 1140 | I8(4,4) = 1.0; |
||
| 1141 | I8(5,5) = 1.0; |
||
| 1142 | I8(6,6) = 1.0; |
||
| 1143 | I8(7,7) = 1.0; |
||
| 1144 | |||
| 1145 | G = I8 - dNmod*mNodeCrd; |
||
| 1146 | |||
| 1147 | // compute gamma vectors |
||
| 1148 | gst = G*hst; |
||
| 1149 | gut = G*hut; |
||
| 1150 | gus = G*hus; |
||
| 1151 | gstu = G*hstu; |
||
| 1152 | |||
| 1153 | // define kinematic and mapping matrices |
||
| 1154 | Bnot.Zero(); |
||
| 1155 | Mben.Zero(); |
||
| 1156 | for (int i = 0; i < 8; i++) { |
||
| 1157 | Bnot(0,3*i) = dNmod(i,0); |
||
| 1158 | Bnot(1,3*i+1) = dNmod(i,1); |
||
| 1159 | Bnot(2,3*i+2) = dNmod(i,2); |
||
| 1160 | Bnot(3,3*i) = dNmod(i,1); |
||
| 1161 | Bnot(3,3*i+1) = dNmod(i,0); |
||
| 1162 | Bnot(4,3*i+1) = dNmod(i,2); |
||
| 1163 | Bnot(4,3*i+2) = dNmod(i,1); |
||
| 1164 | Bnot(5,3*i) = dNmod(i,2); |
||
| 1165 | Bnot(5,3*i+2) = dNmod(i,0); |
||
| 1166 | |||
| 1167 | Mben(0,3*i) = gst(i); |
||
| 1168 | Mben(1,3*i+1) = gst(i); |
||
| 1169 | Mben(2,3*i+2) = gst(i); |
||
| 1170 | Mben(3,3*i) = gut(i); |
||
| 1171 | Mben(4,3*i+1) = gut(i); |
||
| 1172 | Mben(5,3*i+2) = gut(i); |
||
| 1173 | Mben(6,3*i) = gus(i); |
||
| 1174 | Mben(7,3*i+1) = gus(i); |
||
| 1175 | Mben(8,3*i+2) = gus(i); |
||
| 1176 | Mben(9,3*i) = gstu(i); |
||
| 1177 | Mben(10,3*i+1) = gstu(i); |
||
| 1178 | Mben(11,3*i+2) = gstu(i); |
||
| 1179 | } |
||
| 1180 | |||
| 1181 | // define terms for FCF matrix |
||
| 1182 | double HstXX = J0879*Jinv(0,0)*Jinv(0,0) + J0789*Jinv(1,0)*Jinv(1,0) + J619*(Jinv(0,0)*Jinv(1,0) + Jinv(1,0)*Jinv(0,0)); |
||
| 1183 | double HstXY = J0879*Jinv(0,0)*Jinv(0,1) + J0789*Jinv(1,0)*Jinv(1,1) + J619*(Jinv(0,0)*Jinv(1,1) + Jinv(1,0)*Jinv(0,1)); |
||
| 1184 | double HstXZ = J0879*Jinv(0,0)*Jinv(0,2) + J0789*Jinv(1,0)*Jinv(1,2) + J619*(Jinv(0,0)*Jinv(1,2) + Jinv(1,0)*Jinv(0,2)); |
||
| 1185 | double HstYY = J0879*Jinv(0,1)*Jinv(0,1) + J0789*Jinv(1,1)*Jinv(1,1) + J619*(Jinv(0,1)*Jinv(1,1) + Jinv(1,1)*Jinv(0,1)); |
||
| 1186 | double HstYZ = J0879*Jinv(0,1)*Jinv(0,2) + J0789*Jinv(1,1)*Jinv(1,2) + J619*(Jinv(0,1)*Jinv(1,2) + Jinv(1,1)*Jinv(0,2)); |
||
| 1187 | double HstZZ = J0879*Jinv(0,2)*Jinv(0,2) + J0789*Jinv(1,2)*Jinv(1,2) + J619*(Jinv(0,2)*Jinv(1,2) + Jinv(1,2)*Jinv(0,2)); |
||
| 1188 | |||
| 1189 | double HutXX = J0978*Jinv(1,0)*Jinv(1,0) + J0879*Jinv(2,0)*Jinv(2,0) + J417*(Jinv(1,0)*Jinv(2,0) + Jinv(2,0)*Jinv(1,0)); |
||
| 1190 | double HutXY = J0978*Jinv(1,0)*Jinv(1,1) + J0879*Jinv(2,0)*Jinv(2,1) + J417*(Jinv(1,0)*Jinv(2,1) + Jinv(2,0)*Jinv(1,1)); |
||
| 1191 | double HutXZ = J0978*Jinv(1,0)*Jinv(1,2) + J0879*Jinv(2,0)*Jinv(2,2) + J417*(Jinv(1,0)*Jinv(2,2) + Jinv(2,0)*Jinv(1,2)); |
||
| 1192 | double HutYY = J0978*Jinv(1,1)*Jinv(1,1) + J0879*Jinv(2,1)*Jinv(2,1) + J417*(Jinv(1,1)*Jinv(2,1) + Jinv(2,1)*Jinv(1,1)); |
||
| 1193 | double HutYZ = J0978*Jinv(1,1)*Jinv(1,2) + J0879*Jinv(2,1)*Jinv(2,2) + J417*(Jinv(1,1)*Jinv(2,2) + Jinv(2,1)*Jinv(1,2)); |
||
| 1194 | double HutZZ = J0978*Jinv(1,2)*Jinv(1,2) + J0879*Jinv(2,2)*Jinv(2,2) + J417*(Jinv(1,2)*Jinv(2,2) + Jinv(2,2)*Jinv(1,2)); |
||
| 1195 | |||
| 4925 | parduino | 1196 | double HusXX = J0978*Jinv(0,0)*Jinv(0,0) + J0789*Jinv(2,0)*Jinv(2,0) + J518*(Jinv(0,0)*Jinv(2,0) + Jinv(2,0)*Jinv(0,0)); |
| 1197 | double HusXY = J0978*Jinv(0,0)*Jinv(0,1) + J0789*Jinv(2,0)*Jinv(2,1) + J518*(Jinv(0,0)*Jinv(2,1) + Jinv(2,0)*Jinv(0,1)); |
||
| 1198 | double HusXZ = J0978*Jinv(0,0)*Jinv(0,2) + J0789*Jinv(2,0)*Jinv(2,2) + J518*(Jinv(0,0)*Jinv(2,2) + Jinv(2,0)*Jinv(0,2)); |
||
| 1199 | double HusYY = J0978*Jinv(0,1)*Jinv(0,1) + J0789*Jinv(2,1)*Jinv(2,1) + J518*(Jinv(0,1)*Jinv(2,1) + Jinv(2,1)*Jinv(0,1)); |
||
| 1200 | double HusYZ = J0978*Jinv(0,1)*Jinv(0,2) + J0789*Jinv(2,1)*Jinv(2,2) + J518*(Jinv(0,1)*Jinv(2,2) + Jinv(2,1)*Jinv(0,2)); |
||
| 1201 | double HusZZ = J0978*Jinv(0,2)*Jinv(0,2) + J0789*Jinv(2,2)*Jinv(2,2) + J518*(Jinv(0,2)*Jinv(2,2) + Jinv(2,2)*Jinv(0,2)); |
||
| 4654 | parduino | 1202 | |
| 1203 | double HstuXX = J897*Jinv(0,0)*Jinv(0,0) + J798*Jinv(1,0)*Jinv(1,0) + J789*Jinv(2,0)*Jinv(2,0) + J185*(Jinv(0,0)*Jinv(2,0) + Jinv(2,0)*Jinv(0,0)) |
||
| 1204 | + J196*(Jinv(0,0)*Jinv(1,0) + Jinv(1,0)*Jinv(0,0)) + J174*(Jinv(1,0)*Jinv(2,0) + Jinv(2,0)*Jinv(1,0)); |
||
| 1205 | double HstuXY = J897*Jinv(0,0)*Jinv(0,1) + J798*Jinv(1,0)*Jinv(1,1) + J789*Jinv(2,0)*Jinv(2,1) + J185*(Jinv(0,0)*Jinv(2,1) + Jinv(2,0)*Jinv(0,1)) |
||
| 1206 | + J196*(Jinv(0,0)*Jinv(1,1) + Jinv(1,0)*Jinv(0,1)) + J174*(Jinv(1,0)*Jinv(2,1) + Jinv(2,0)*Jinv(1,1)); |
||
| 1207 | double HstuXZ = J897*Jinv(0,0)*Jinv(0,2) + J798*Jinv(1,0)*Jinv(1,2) + J789*Jinv(2,0)*Jinv(2,2) + J185*(Jinv(0,0)*Jinv(2,2) + Jinv(2,0)*Jinv(0,2)) |
||
| 1208 | + J196*(Jinv(0,0)*Jinv(1,2) + Jinv(1,0)*Jinv(0,2)) + J174*(Jinv(1,0)*Jinv(2,2) + Jinv(2,0)*Jinv(1,2)); |
||
| 1209 | double HstuYY = J897*Jinv(0,1)*Jinv(0,1) + J798*Jinv(1,1)*Jinv(1,1) + J789*Jinv(2,1)*Jinv(2,1) + J185*(Jinv(0,1)*Jinv(2,1) + Jinv(2,1)*Jinv(0,1)) |
||
| 1210 | + J196*(Jinv(0,1)*Jinv(1,1) + Jinv(1,1)*Jinv(0,1)) + J174*(Jinv(1,1)*Jinv(2,1) + Jinv(2,1)*Jinv(1,1)); |
||
| 1211 | double HstuYZ = J897*Jinv(0,1)*Jinv(0,2) + J798*Jinv(1,1)*Jinv(1,2) + J789*Jinv(2,1)*Jinv(2,2) + J185*(Jinv(0,1)*Jinv(2,2) + Jinv(2,1)*Jinv(0,2)) |
||
| 1212 | + J196*(Jinv(0,1)*Jinv(1,2) + Jinv(1,1)*Jinv(0,2)) + J174*(Jinv(1,1)*Jinv(2,2) + Jinv(2,1)*Jinv(1,2)); |
||
| 1213 | double HstuZZ = J897*Jinv(0,2)*Jinv(0,2) + J798*Jinv(1,2)*Jinv(1,2) + J789*Jinv(2,2)*Jinv(2,2) + J185*(Jinv(0,2)*Jinv(2,2) + Jinv(2,2)*Jinv(0,2)) |
||
| 1214 | + J196*(Jinv(0,2)*Jinv(1,2) + Jinv(1,2)*Jinv(0,2)) + J174*(Jinv(1,2)*Jinv(2,2) + Jinv(2,2)*Jinv(1,2)); |
||
| 1215 | |||
| 1216 | double IttXX = J0879*Jinv(0,0)*Jinv(2,0) + J417*Jinv(0,0)*Jinv(1,0) + J518*Jinv(1,0)*Jinv(1,0) + J619*Jinv(1,0)*Jinv(2,0); |
||
| 1217 | double IttXY = J0879*Jinv(0,0)*Jinv(2,1) + J417*Jinv(0,0)*Jinv(1,1) + J518*Jinv(1,0)*Jinv(1,1) + J619*Jinv(1,0)*Jinv(2,1); |
||
| 1218 | double IttXZ = J0879*Jinv(0,0)*Jinv(2,2) + J417*Jinv(0,0)*Jinv(1,2) + J518*Jinv(1,0)*Jinv(1,2) + J619*Jinv(1,0)*Jinv(2,2); |
||
| 1219 | double IttYX = J0879*Jinv(0,1)*Jinv(2,0) + J417*Jinv(0,1)*Jinv(1,0) + J518*Jinv(1,1)*Jinv(1,0) + J619*Jinv(1,1)*Jinv(2,0); |
||
| 1220 | double IttYY = J0879*Jinv(0,1)*Jinv(2,1) + J417*Jinv(0,1)*Jinv(1,1) + J518*Jinv(1,1)*Jinv(1,1) + J619*Jinv(1,1)*Jinv(2,1); |
||
| 1221 | double IttYZ = J0879*Jinv(0,1)*Jinv(2,2) + J417*Jinv(0,1)*Jinv(1,2) + J518*Jinv(1,1)*Jinv(1,2) + J619*Jinv(1,1)*Jinv(2,2); |
||
| 1222 | double IttZX = J0879*Jinv(0,2)*Jinv(2,0) + J417*Jinv(0,2)*Jinv(1,0) + J518*Jinv(1,2)*Jinv(1,0) + J619*Jinv(1,2)*Jinv(2,0); |
||
| 1223 | double IttZY = J0879*Jinv(0,2)*Jinv(2,1) + J417*Jinv(0,2)*Jinv(1,1) + J518*Jinv(1,2)*Jinv(1,1) + J619*Jinv(1,2)*Jinv(2,1); |
||
| 1224 | double IttZZ = J0879*Jinv(0,2)*Jinv(2,2) + J417*Jinv(0,2)*Jinv(1,2) + J518*Jinv(1,2)*Jinv(1,2) + J619*Jinv(1,2)*Jinv(2,2); |
||
| 1225 | |||
| 4925 | parduino | 1226 | double IssXX = J0789*Jinv(1,0)*Jinv(2,0) + J417*Jinv(0,0)*Jinv(0,0) + J518*Jinv(1,0)*Jinv(0,0) + J619*Jinv(0,0)*Jinv(2,0); |
| 4654 | parduino | 1227 | double IssXY = J0789*Jinv(1,0)*Jinv(2,1) + J417*Jinv(0,0)*Jinv(0,1) + J518*Jinv(1,0)*Jinv(0,1) + J619*Jinv(0,0)*Jinv(2,1); |
| 1228 | double IssXZ = J0789*Jinv(1,0)*Jinv(2,2) + J417*Jinv(0,0)*Jinv(0,2) + J518*Jinv(1,0)*Jinv(0,2) + J619*Jinv(0,0)*Jinv(2,2); |
||
| 1229 | double IssYX = J0789*Jinv(1,1)*Jinv(2,0) + J417*Jinv(0,1)*Jinv(0,0) + J518*Jinv(1,1)*Jinv(0,0) + J619*Jinv(0,1)*Jinv(2,0); |
||
| 1230 | double IssYY = J0789*Jinv(1,1)*Jinv(2,1) + J417*Jinv(0,1)*Jinv(0,1) + J518*Jinv(1,1)*Jinv(0,1) + J619*Jinv(0,1)*Jinv(2,1); |
||
| 1231 | double IssYZ = J0789*Jinv(1,1)*Jinv(2,2) + J417*Jinv(0,1)*Jinv(0,2) + J518*Jinv(1,1)*Jinv(0,2) + J619*Jinv(0,1)*Jinv(2,2); |
||
| 1232 | double IssZX = J0789*Jinv(1,2)*Jinv(2,0) + J417*Jinv(0,2)*Jinv(0,0) + J518*Jinv(1,2)*Jinv(0,0) + J619*Jinv(0,2)*Jinv(2,0); |
||
| 1233 | double IssZY = J0789*Jinv(1,2)*Jinv(2,1) + J417*Jinv(0,2)*Jinv(0,1) + J518*Jinv(1,2)*Jinv(0,1) + J619*Jinv(0,2)*Jinv(2,1); |
||
| 1234 | double IssZZ = J0789*Jinv(1,2)*Jinv(2,2) + J417*Jinv(0,2)*Jinv(0,2) + J518*Jinv(1,2)*Jinv(0,2) + J619*Jinv(0,2)*Jinv(2,2); |
||
| 1235 | |||
| 1236 | double IuuXX = J0978*Jinv(1,0)*Jinv(0,0) + J417*Jinv(2,0)*Jinv(0,0) + J518*Jinv(1,0)*Jinv(2,0) + J619*Jinv(2,0)*Jinv(2,0); |
||
| 1237 | double IuuXY = J0978*Jinv(1,0)*Jinv(0,1) + J417*Jinv(2,0)*Jinv(0,1) + J518*Jinv(1,0)*Jinv(2,1) + J619*Jinv(2,0)*Jinv(2,1); |
||
| 1238 | double IuuXZ = J0978*Jinv(1,0)*Jinv(0,2) + J417*Jinv(2,0)*Jinv(0,2) + J518*Jinv(1,0)*Jinv(2,2) + J619*Jinv(2,0)*Jinv(2,2); |
||
| 1239 | double IuuYX = J0978*Jinv(1,1)*Jinv(0,0) + J417*Jinv(2,1)*Jinv(0,0) + J518*Jinv(1,1)*Jinv(2,0) + J619*Jinv(2,1)*Jinv(2,0); |
||
| 1240 | double IuuYY = J0978*Jinv(1,1)*Jinv(0,1) + J417*Jinv(2,1)*Jinv(0,1) + J518*Jinv(1,1)*Jinv(2,1) + J619*Jinv(2,1)*Jinv(2,1); |
||
| 1241 | double IuuYZ = J0978*Jinv(1,1)*Jinv(0,2) + J417*Jinv(2,1)*Jinv(0,2) + J518*Jinv(1,1)*Jinv(2,2) + J619*Jinv(2,1)*Jinv(2,2); |
||
| 1242 | double IuuZX = J0978*Jinv(1,2)*Jinv(0,0) + J417*Jinv(2,2)*Jinv(0,0) + J518*Jinv(1,2)*Jinv(2,0) + J619*Jinv(2,2)*Jinv(2,0); |
||
| 1243 | double IuuZY = J0978*Jinv(1,2)*Jinv(0,1) + J417*Jinv(2,2)*Jinv(0,1) + J518*Jinv(1,2)*Jinv(2,1) + J619*Jinv(2,2)*Jinv(2,1); |
||
| 1244 | double IuuZZ = J0978*Jinv(1,2)*Jinv(0,2) + J417*Jinv(2,2)*Jinv(0,2) + J518*Jinv(1,2)*Jinv(2,2) + J619*Jinv(2,2)*Jinv(2,2); |
||
| 4925 | parduino | 1245 | |
| 1246 | double IstXX = J31013*Jinv(0,0)*Jinv(0,0) + J31310*Jinv(1,0)*Jinv(1,0) + J21411*Jinv(2,0)*Jinv(1,0) + J11512*Jinv(2,0)*Jinv(0,0) + J16*(Jinv(0,0)*Jinv(1,0) + Jinv(1,0)*Jinv(0,0)); |
||
| 1247 | double IstXY = J31013*Jinv(0,0)*Jinv(0,1) + J31310*Jinv(1,0)*Jinv(1,1) + J21411*Jinv(2,0)*Jinv(1,1) + J11512*Jinv(2,0)*Jinv(0,1) + J16*(Jinv(0,0)*Jinv(1,1) + Jinv(1,0)*Jinv(0,1)); |
||
| 1248 | double IstXZ = J31013*Jinv(0,0)*Jinv(0,2) + J31310*Jinv(1,0)*Jinv(1,2) + J21411*Jinv(2,0)*Jinv(1,2) + J11512*Jinv(2,0)*Jinv(0,2) + J16*(Jinv(0,0)*Jinv(1,2) + Jinv(1,0)*Jinv(0,2)); |
||
| 1249 | double IstYX = J31013*Jinv(0,1)*Jinv(0,0) + J31310*Jinv(1,1)*Jinv(1,0) + J21411*Jinv(2,1)*Jinv(1,0) + J11512*Jinv(2,1)*Jinv(0,0) + J16*(Jinv(0,1)*Jinv(1,0) + Jinv(1,1)*Jinv(0,0)); |
||
| 1250 | double IstYY = J31013*Jinv(0,1)*Jinv(0,1) + J31310*Jinv(1,1)*Jinv(1,1) + J21411*Jinv(2,1)*Jinv(1,1) + J11512*Jinv(2,1)*Jinv(0,1) + J16*(Jinv(0,1)*Jinv(1,1) + Jinv(1,1)*Jinv(0,1)); |
||
| 1251 | double IstYZ = J31013*Jinv(0,1)*Jinv(0,2) + J31310*Jinv(1,1)*Jinv(1,2) + J21411*Jinv(2,1)*Jinv(1,2) + J11512*Jinv(2,1)*Jinv(0,2) + J16*(Jinv(0,1)*Jinv(1,2) + Jinv(1,1)*Jinv(0,2)); |
||
| 1252 | double IstZX = J31013*Jinv(0,2)*Jinv(0,0) + J31310*Jinv(1,2)*Jinv(1,0) + J21411*Jinv(2,2)*Jinv(1,0) + J11512*Jinv(2,2)*Jinv(0,0) + J16*(Jinv(0,2)*Jinv(1,0) + Jinv(1,2)*Jinv(0,0)); |
||
| 1253 | double IstZY = J31013*Jinv(0,2)*Jinv(0,1) + J31310*Jinv(1,2)*Jinv(1,1) + J21411*Jinv(2,2)*Jinv(1,1) + J11512*Jinv(2,2)*Jinv(0,1) + J16*(Jinv(0,2)*Jinv(1,1) + Jinv(1,2)*Jinv(0,1)); |
||
| 1254 | double IstZZ = J31013*Jinv(0,2)*Jinv(0,2) + J31310*Jinv(1,2)*Jinv(1,2) + J21411*Jinv(2,2)*Jinv(1,2) + J11512*Jinv(2,2)*Jinv(0,2) + J16*(Jinv(0,2)*Jinv(1,2) + Jinv(1,2)*Jinv(0,2)); |
||
| 4654 | parduino | 1255 | |
| 4925 | parduino | 1256 | double IutXX = J21114*Jinv(0,0)*Jinv(1,0) + J31013*Jinv(0,0)*Jinv(2,0) + J11215*Jinv(1,0)*Jinv(1,0) + J11512*Jinv(2,0)*Jinv(2,0) + J16*(Jinv(1,0)*Jinv(2,0) + Jinv(2,0)*Jinv(1,0)); |
| 1257 | double IutXY = J21114*Jinv(0,0)*Jinv(1,1) + J31013*Jinv(0,0)*Jinv(2,1) + J11215*Jinv(1,0)*Jinv(1,1) + J11512*Jinv(2,0)*Jinv(2,1) + J16*(Jinv(1,0)*Jinv(2,1) + Jinv(2,0)*Jinv(1,1)); |
||
| 1258 | double IutXZ = J21114*Jinv(0,0)*Jinv(1,2) + J31013*Jinv(0,0)*Jinv(2,2) + J11215*Jinv(1,0)*Jinv(1,2) + J11512*Jinv(2,0)*Jinv(2,2) + J16*(Jinv(1,0)*Jinv(2,2) + Jinv(2,0)*Jinv(1,2)); |
||
| 1259 | double IutYX = J21114*Jinv(0,1)*Jinv(1,0) + J31013*Jinv(0,1)*Jinv(2,0) + J11215*Jinv(1,1)*Jinv(1,0) + J11512*Jinv(2,1)*Jinv(2,0) + J16*(Jinv(1,1)*Jinv(2,0) + Jinv(2,1)*Jinv(1,0)); |
||
| 1260 | double IutYY = J21114*Jinv(0,1)*Jinv(1,1) + J31013*Jinv(0,1)*Jinv(2,1) + J11215*Jinv(1,1)*Jinv(1,1) + J11512*Jinv(2,1)*Jinv(2,1) + J16*(Jinv(1,1)*Jinv(2,1) + Jinv(2,1)*Jinv(1,1)); |
||
| 1261 | double IutYZ = J21114*Jinv(0,1)*Jinv(1,2) + J31013*Jinv(0,1)*Jinv(2,2) + J11215*Jinv(1,1)*Jinv(1,2) + J11512*Jinv(2,1)*Jinv(2,2) + J16*(Jinv(1,1)*Jinv(2,2) + Jinv(2,1)*Jinv(1,2)); |
||
| 1262 | double IutZX = J21114*Jinv(0,2)*Jinv(1,0) + J31013*Jinv(0,2)*Jinv(2,0) + J11215*Jinv(1,2)*Jinv(1,0) + J11512*Jinv(2,2)*Jinv(2,0) + J16*(Jinv(1,2)*Jinv(2,0) + Jinv(2,2)*Jinv(1,0)); |
||
| 1263 | double IutZY = J21114*Jinv(0,2)*Jinv(1,1) + J31013*Jinv(0,2)*Jinv(2,1) + J11215*Jinv(1,2)*Jinv(1,1) + J11512*Jinv(2,2)*Jinv(2,1) + J16*(Jinv(1,2)*Jinv(2,1) + Jinv(2,2)*Jinv(1,1)); |
||
| 1264 | double IutZZ = J21114*Jinv(0,2)*Jinv(1,2) + J31013*Jinv(0,2)*Jinv(2,2) + J11215*Jinv(1,2)*Jinv(1,2) + J11512*Jinv(2,2)*Jinv(2,2) + J16*(Jinv(1,2)*Jinv(2,2) + Jinv(2,2)*Jinv(1,2)); |
||
| 4654 | parduino | 1265 | |
| 4925 | parduino | 1266 | double IusXX = J21114*Jinv(0,0)*Jinv(0,0) + J11215*Jinv(1,0)*Jinv(0,0) + J31310*Jinv(1,0)*Jinv(2,0) + J21411*Jinv(2,0)*Jinv(2,0) + J16*(Jinv(0,0)*Jinv(2,0) + Jinv(2,0)*Jinv(0,0)); |
| 1267 | double IusXY = J21114*Jinv(0,0)*Jinv(0,1) + J11215*Jinv(1,0)*Jinv(0,1) + J31310*Jinv(1,0)*Jinv(2,1) + J21411*Jinv(2,0)*Jinv(2,1) + J16*(Jinv(0,0)*Jinv(2,1) + Jinv(2,0)*Jinv(0,1)); |
||
| 1268 | double IusXZ = J21114*Jinv(0,0)*Jinv(0,2) + J11215*Jinv(1,0)*Jinv(0,2) + J31310*Jinv(1,0)*Jinv(2,2) + J21411*Jinv(2,0)*Jinv(2,2) + J16*(Jinv(0,0)*Jinv(2,2) + Jinv(2,0)*Jinv(0,2)); |
||
| 1269 | double IusYX = J21114*Jinv(0,1)*Jinv(0,0) + J11215*Jinv(1,1)*Jinv(0,0) + J31310*Jinv(1,1)*Jinv(2,0) + J21411*Jinv(2,1)*Jinv(2,0) + J16*(Jinv(0,1)*Jinv(2,0) + Jinv(2,1)*Jinv(0,0)); |
||
| 1270 | double IusYY = J21114*Jinv(0,1)*Jinv(0,1) + J11215*Jinv(1,1)*Jinv(0,1) + J31310*Jinv(1,1)*Jinv(2,1) + J21411*Jinv(2,1)*Jinv(2,1) + J16*(Jinv(0,1)*Jinv(2,1) + Jinv(2,1)*Jinv(0,1)); |
||
| 1271 | double IusYZ = J21114*Jinv(0,1)*Jinv(0,2) + J11215*Jinv(1,1)*Jinv(0,2) + J31310*Jinv(1,1)*Jinv(2,2) + J21411*Jinv(2,1)*Jinv(2,2) + J16*(Jinv(0,1)*Jinv(2,2) + Jinv(2,1)*Jinv(0,2)); |
||
| 1272 | double IusZX = J21114*Jinv(0,2)*Jinv(0,0) + J11215*Jinv(1,2)*Jinv(0,0) + J31310*Jinv(1,2)*Jinv(2,0) + J21411*Jinv(2,2)*Jinv(2,0) + J16*(Jinv(0,2)*Jinv(2,0) + Jinv(2,2)*Jinv(0,0)); |
||
| 1273 | double IusZY = J21114*Jinv(0,2)*Jinv(0,1) + J11215*Jinv(1,2)*Jinv(0,1) + J31310*Jinv(1,2)*Jinv(2,1) + J21411*Jinv(2,2)*Jinv(2,1) + J16*(Jinv(0,2)*Jinv(2,1) + Jinv(2,2)*Jinv(0,1)); |
||
| 1274 | double IusZZ = J21114*Jinv(0,2)*Jinv(0,2) + J11215*Jinv(1,2)*Jinv(0,2) + J31310*Jinv(1,2)*Jinv(2,2) + J21411*Jinv(2,2)*Jinv(2,2) + J16*(Jinv(0,2)*Jinv(2,2) + Jinv(2,2)*Jinv(0,2)); |
||
| 4654 | parduino | 1275 | |
| 4925 | parduino | 1276 | // constitutive constants from material |
| 4654 | parduino | 1277 | const Matrix &CmatI = theMaterial->getInitialTangent(); |
| 4656 | parduino | 1278 | double C1 = CmatI(0,0); |
| 1279 | double C2 = CmatI(0,1); |
||
| 1280 | double C3 = CmatI(3,3); |
||
| 4925 | parduino | 1281 | double C4 = C2 + C3; |
| 4654 | parduino | 1282 | |
| 1283 | // define the 12x12 matrix FCF in 3x3 blocks |
||
| 4656 | parduino | 1284 | // block11 |
| 4654 | parduino | 1285 | FCF(0,0) = C1*HstXX + C3*(HstYY + HstZZ); |
| 4925 | parduino | 1286 | FCF(0,1) = C4*HstXY; |
| 1287 | FCF(0,2) = C4*HstXZ; |
||
| 1288 | FCF(1,0) = FCF(0,1); |
||
| 4654 | parduino | 1289 | FCF(1,1) = C1*HstYY + C3*(HstXX + HstZZ); |
| 4925 | parduino | 1290 | FCF(1,2) = C4*HstYZ; |
| 1291 | FCF(2,0) = FCF(0,2); |
||
| 1292 | FCF(2,1) = FCF(1,2); |
||
| 4654 | parduino | 1293 | FCF(2,2) = C1*HstZZ + C3*(HstYY + HstXX); |
| 1294 | |||
| 1295 | // block12 |
||
| 1296 | FCF(0,3) = C1*IttXX + C3*(IttYY + IttZZ); |
||
| 1297 | FCF(0,4) = C2*IttXY + C3*IttYX; |
||
| 1298 | FCF(0,5) = C2*IttXZ + C3*IttZX; |
||
| 1299 | FCF(1,3) = C2*IttYX + C3*IttXY; |
||
| 1300 | FCF(1,4) = C1*IttYY + C3*(IttXX + IttZZ); |
||
| 1301 | FCF(1,5) = C2*IttYZ + C3*IttZY; |
||
| 1302 | FCF(2,3) = C2*IttZX + C3*IttXZ; |
||
| 1303 | FCF(2,4) = C2*IttZY + C3*IttYZ; |
||
| 1304 | FCF(2,5) = C1*IttZZ + C3*(IttYY + IttXX); |
||
| 1305 | |||
| 1306 | // block13 |
||
| 1307 | FCF(0,6) = C1*IssXX + C3*(IssYY + IssZZ); |
||
| 1308 | FCF(0,7) = C2*IssXY + C3*IssYX; |
||
| 1309 | FCF(0,8) = C2*IssXZ + C3*IssZX; |
||
| 1310 | FCF(1,6) = C2*IssYX + C3*IssXY; |
||
| 1311 | FCF(1,7) = C1*IssYY + C3*(IssXX + IssZZ); |
||
| 1312 | FCF(1,8) = C2*IssYZ + C3*IssZY; |
||
| 1313 | FCF(2,6) = C2*IssZX + C3*IssXZ; |
||
| 1314 | FCF(2,7) = C2*IssZY + C3*IssYZ; |
||
| 1315 | FCF(2,8) = C1*IssZZ + C3*(IssYY + IssXX); |
||
| 1316 | |||
| 4659 | parduino | 1317 | /*// block14 |
| 4654 | parduino | 1318 | FCF(0,9) = C1*IstXX + C3*(IstYY + IstZZ); |
| 1319 | FCF(0,10) = C2*IstYX + C3*IstXY; |
||
| 1320 | FCF(0,11) = C2*IstZX + C3*IstXZ; |
||
| 1321 | FCF(1,9) = C2*IstXY + C3*IstYX; |
||
| 1322 | FCF(1,10) = C1*IstYY + C3*(IstXX + IstZZ); |
||
| 1323 | FCF(1,11) = C2*IstZY + C3*IstYZ; |
||
| 1324 | FCF(2,9) = C2*IstXZ + C3*IstZX; |
||
| 1325 | FCF(2,10) = C2*IstYZ + C3*IstZY; |
||
| 4659 | parduino | 1326 | FCF(2,11) = C1*IstZZ + C3*(IstYY + IstXX);*/ |
| 4654 | parduino | 1327 | |
| 4659 | parduino | 1328 | // block14 |
| 1329 | FCF(0,9) = C3*(IstYY + IstZZ); |
||
| 1330 | FCF(0,10) = C3*IstXY; |
||
| 1331 | FCF(0,11) = C3*IstXZ; |
||
| 1332 | FCF(1,9) = C3*IstYX; |
||
| 1333 | FCF(1,10) = C3*(IstXX + IstZZ); |
||
| 1334 | FCF(1,11) = C3*IstYZ; |
||
| 1335 | FCF(2,9) = C3*IstZX; |
||
| 1336 | FCF(2,10) = C3*IstZY; |
||
| 1337 | FCF(2,11) = C3*(IstYY + IstXX); |
||
| 1338 | |||
| 4654 | parduino | 1339 | // block21 |
| 1340 | FCF(3,0) = C1*IttXX + C3*(IttYY + IttZZ); |
||
| 1341 | FCF(3,1) = C2*IttYX + C3*IttXY; |
||
| 1342 | FCF(3,2) = C2*IttZX + C3*IttXZ; |
||
| 1343 | FCF(4,0) = C2*IttXY + C3*IttYX; |
||
| 1344 | FCF(4,1) = C1*IttYY + C3*(IttXX + IttZZ); |
||
| 1345 | FCF(4,2) = C2*IttZY + C3*IttYZ; |
||
| 1346 | FCF(5,0) = C2*IttXZ + C3*IttZX; |
||
| 1347 | FCF(5,1) = C2*IttYZ + C3*IttZY; |
||
| 1348 | FCF(5,2) = C1*IttZZ + C3*(IttYY + IttXX); |
||
| 1349 | |||
| 1350 | // block22 |
||
| 1351 | FCF(3,3) = C1*HutXX + C3*(HutYY + HutZZ); |
||
| 4925 | parduino | 1352 | FCF(3,4) = C4*HutXY; |
| 1353 | FCF(3,5) = C4*HutXZ; |
||
| 1354 | FCF(4,3) = FCF(3,4); |
||
| 4654 | parduino | 1355 | FCF(4,4) = C1*HutYY + C3*(HutXX + HutZZ); |
| 4925 | parduino | 1356 | FCF(4,5) = C4*HutYZ; |
| 1357 | FCF(5,3) = FCF(3,5); |
||
| 1358 | FCF(5,4) = FCF(4,5); |
||
| 4654 | parduino | 1359 | FCF(5,5) = C1*HutZZ + C3*(HutYY + HutXX); |
| 1360 | |||
| 1361 | // block23 |
||
| 1362 | FCF(3,6) = C1*IuuXX + C3*(IuuYY + IuuZZ); |
||
| 1363 | FCF(3,7) = C2*IuuXY + C3*IuuYX; |
||
| 1364 | FCF(3,8) = C2*IuuXZ + C3*IuuZX; |
||
| 1365 | FCF(4,6) = C2*IuuYX + C3*IuuXY; |
||
| 1366 | FCF(4,7) = C1*IuuYY + C3*(IuuXX + IuuZZ); |
||
| 1367 | FCF(4,8) = C2*IuuYZ + C3*IuuZY; |
||
| 1368 | FCF(5,6) = C2*IuuZX + C3*IuuXZ; |
||
| 1369 | FCF(5,7) = C2*IuuZY + C3*IuuYZ; |
||
| 1370 | FCF(5,8) = C1*IuuZZ + C3*(IuuYY + IuuXX); |
||
| 1371 | |||
| 4659 | parduino | 1372 | /*// block24 |
| 4654 | parduino | 1373 | FCF(3,9) = C1*IutXX + C3*(IutYY + IutZZ); |
| 1374 | FCF(3,10) = C2*IutYX + C3*IutXY; |
||
| 1375 | FCF(3,11) = C2*IutZX + C3*IutXZ; |
||
| 1376 | FCF(4,9) = C2*IutXY + C3*IutYX; |
||
| 1377 | FCF(4,10) = C1*IutYY + C3*(IutXX + IutZZ); |
||
| 1378 | FCF(4,11) = C2*IutZY + C3*IutYZ; |
||
| 1379 | FCF(5,9) = C2*IutXZ + C3*IutZX; |
||
| 1380 | FCF(5,10) = C2*IutYZ + C3*IutZY; |
||
| 4659 | parduino | 1381 | FCF(5,11) = C1*IutZZ + C3*(IutYY + IutXX);*/ |
| 4654 | parduino | 1382 | |
| 4659 | parduino | 1383 | // block24 |
| 1384 | FCF(3,9) = C3*(IutYY + IutZZ); |
||
| 1385 | FCF(3,10) = C3*IutXY; |
||
| 1386 | FCF(3,11) = C3*IutXZ; |
||
| 1387 | FCF(4,9) = C3*IutYX; |
||
| 1388 | FCF(4,10) = C3*(IutXX + IutZZ); |
||
| 1389 | FCF(4,11) = C3*IutYZ; |
||
| 1390 | FCF(5,9) = C3*IutZX; |
||
| 1391 | FCF(5,10) = C3*IutZY; |
||
| 1392 | FCF(5,11) = C3*(IutYY + IutXX); |
||
| 1393 | |||
| 4654 | parduino | 1394 | // block31 |
| 1395 | FCF(6,0) = C1*IssXX + C3*(IssYY + IssZZ); |
||
| 1396 | FCF(6,1) = C2*IssYX + C3*IssXY; |
||
| 1397 | FCF(6,2) = C2*IssZX + C3*IssXZ; |
||
| 1398 | FCF(7,0) = C2*IssXY + C3*IssYX; |
||
| 1399 | FCF(7,1) = C1*IssYY + C3*(IssXX + IssZZ); |
||
| 1400 | FCF(7,2) = C2*IssZY + C3*IssYZ; |
||
| 1401 | FCF(8,0) = C2*IssXZ + C3*IssZX; |
||
| 1402 | FCF(8,1) = C2*IssYZ + C3*IssZY; |
||
| 1403 | FCF(8,2) = C1*IssZZ + C3*(IssYY + IssXX); |
||
| 1404 | |||
| 1405 | // block32 |
||
| 1406 | FCF(6,3) = C1*IuuXX + C3*(IuuYY + IuuZZ); |
||
| 1407 | FCF(6,4) = C2*IuuYX + C3*IuuXY; |
||
| 1408 | FCF(6,5) = C2*IuuZX + C3*IuuXZ; |
||
| 1409 | FCF(7,3) = C2*IuuXY + C3*IuuYX; |
||
| 1410 | FCF(7,4) = C1*IuuYY + C3*(IuuXX + IuuZZ); |
||
| 1411 | FCF(7,5) = C2*IuuZY + C3*IuuYZ; |
||
| 1412 | FCF(8,3) = C2*IuuXZ + C3*IuuZX; |
||
| 1413 | FCF(8,4) = C2*IuuYZ + C3*IuuZY; |
||
| 1414 | FCF(8,5) = C1*IuuZZ + C3*(IuuYY + IuuXX); |
||
| 1415 | |||
| 1416 | // block33 |
||
| 1417 | FCF(6,6) = C1*HusXX + C3*(HusYY + HusZZ); |
||
| 4925 | parduino | 1418 | FCF(6,7) = C4*HusXY; |
| 1419 | FCF(6,8) = C4*HusXZ; |
||
| 1420 | FCF(7,6) = FCF(6,7); |
||
| 4654 | parduino | 1421 | FCF(7,7) = C1*HusYY + C3*(HusXX + HusZZ); |
| 4925 | parduino | 1422 | FCF(7,8) = C4*HusYZ; |
| 1423 | FCF(8,6) = FCF(6,8); |
||
| 1424 | FCF(8,7) = FCF(7,8); |
||
| 4654 | parduino | 1425 | FCF(8,8) = C1*HusZZ + C3*(HusYY + HusXX); |
| 1426 | |||
| 4659 | parduino | 1427 | /*// block34 |
| 4654 | parduino | 1428 | FCF(6,9) = C1*IusXX + C3*(IusYY + IusZZ); |
| 1429 | FCF(6,10) = C2*IusYX + C3*IusXY; |
||
| 1430 | FCF(6,11) = C2*IusZX + C3*IusXZ; |
||
| 1431 | FCF(7,9) = C2*IusXY + C3*IusYX; |
||
| 1432 | FCF(7,10) = C1*IusYY + C3*(IusXX + IusZZ); |
||
| 1433 | FCF(7,11) = C2*IusZY + C3*IusYZ; |
||
| 1434 | FCF(8,9) = C2*IusXZ + C3*IusZX; |
||
| 1435 | FCF(8,10) = C2*IusYZ + C3*IusZY; |
||
| 4659 | parduino | 1436 | FCF(8,11) = C1*IusZZ + C3*(IusYY + IusXX);*/ |
| 4654 | parduino | 1437 | |
| 4659 | parduino | 1438 | // block34 |
| 1439 | FCF(6,9) = C3*(IusYY + IusZZ); |
||
| 1440 | FCF(6,10) = C3*IusXY; |
||
| 1441 | FCF(6,11) = C3*IusXZ; |
||
| 1442 | FCF(7,9) = C3*IusYX; |
||
| 1443 | FCF(7,10) = C3*(IusXX + IusZZ); |
||
| 1444 | FCF(7,11) = C3*IusYZ; |
||
| 1445 | FCF(8,9) = C3*IusZX; |
||
| 1446 | FCF(8,10) = C3*IusZY; |
||
| 1447 | FCF(8,11) = C3*(IusYY + IusXX); |
||
| 1448 | |||
| 1449 | /*// block41 |
||
| 4654 | parduino | 1450 | FCF(9,0) = C1*IstXX + C3*(IstYY + IstZZ); |
| 1451 | FCF(9,1) = C2*IstXY + C3*IstYX; |
||
| 1452 | FCF(9,2) = C2*IstXZ + C3*IstZX; |
||
| 1453 | FCF(10,0) = C2*IstYX + C3*IstXY; |
||
| 1454 | FCF(10,1) = C1*IstYY + C3*(IstXX + IstZZ); |
||
| 1455 | FCF(10,2) = C2*IstYZ + C3*IstZY; |
||
| 1456 | FCF(11,0) = C2*IstZX + C3*IstXZ; |
||
| 1457 | FCF(11,1) = C2*IstZY + C3*IstYZ; |
||
| 1458 | FCF(11,2) = C1*IstZZ + C3*(IstYY + IstXX); |
||
| 1459 | |||
| 1460 | // block42 |
||
| 1461 | FCF(9,3) = C1*IutXX + C3*(IutYY + IutZZ); |
||
| 1462 | FCF(9,4) = C2*IutXY + C3*IutYX; |
||
| 1463 | FCF(9,5) = C2*IutXZ + C3*IutZX; |
||
| 1464 | FCF(10,3) = C2*IutYX + C3*IutXY; |
||
| 1465 | FCF(10,4) = C1*IutYY + C3*(IutXX + IutZZ); |
||
| 1466 | FCF(10,5) = C2*IutYZ + C3*IutZY; |
||
| 1467 | FCF(11,3) = C2*IutZX + C3*IutXZ; |
||
| 1468 | FCF(11,4) = C2*IutZY + C3*IutYZ; |
||
| 1469 | FCF(11,5) = C1*IutZZ + C3*(IutYY + IutXX); |
||
| 1470 | |||
| 1471 | // block43 |
||
| 1472 | FCF(9,6) = C1*IusXX + C3*(IusYY + IusZZ); |
||
| 1473 | FCF(9,7) = C2*IusXY + C3*IusYX; |
||
| 1474 | FCF(9,8) = C2*IusXZ + C3*IusZX; |
||
| 1475 | FCF(10,6) = C2*IusYX + C3*IusXY; |
||
| 1476 | FCF(10,7) = C1*IusYY + C3*(IusXX + IusZZ); |
||
| 1477 | FCF(10,8) = C2*IusYZ + C3*IusZY; |
||
| 1478 | FCF(11,6) = C2*IusZX + C3*IusXZ; |
||
| 1479 | FCF(11,7) = C2*IusZY + C3*IusYZ; |
||
| 4659 | parduino | 1480 | FCF(11,8) = C1*IusZZ + C3*(IusYY + IusXX);*/ |
| 4654 | parduino | 1481 | |
| 4659 | parduino | 1482 | // block41 |
| 1483 | FCF(9,0) = C3*(IstYY + IstZZ); |
||
| 1484 | FCF(9,1) = C3*IstYX; |
||
| 1485 | FCF(9,2) = C3*IstZX; |
||
| 1486 | FCF(10,0) = C3*IstXY; |
||
| 1487 | FCF(10,1) = C3*(IstXX + IstZZ); |
||
| 1488 | FCF(10,2) = C3*IstZY; |
||
| 1489 | FCF(11,0) = C3*IstXZ; |
||
| 1490 | FCF(11,1) = C3*IstYZ; |
||
| 1491 | FCF(11,2) = C3*(IstYY + IstXX); |
||
| 1492 | |||
| 1493 | // block42 |
||
| 1494 | FCF(9,3) = C3*(IutYY + IutZZ); |
||
| 1495 | FCF(9,4) = C3*IutYX; |
||
| 1496 | FCF(9,5) = C3*IutZX; |
||
| 1497 | FCF(10,3) = C3*IutXY; |
||
| 1498 | FCF(10,4) = C3*(IutXX + IutZZ); |
||
| 1499 | FCF(10,5) = C3*IutZY; |
||
| 1500 | FCF(11,3) = C3*IutXZ; |
||
| 1501 | FCF(11,4) = C3*IutYZ; |
||
| 1502 | FCF(11,5) = C3*(IutYY + IutXX); |
||
| 1503 | |||
| 1504 | // block43 |
||
| 1505 | FCF(9,6) = C3*(IusYY + IusZZ); |
||
| 1506 | FCF(9,7) = C3*IusYX; |
||
| 1507 | FCF(9,8) = C3*IusZX; |
||
| 1508 | FCF(10,6) = C3*IusXY; |
||
| 1509 | FCF(10,7) = C3*(IusXX + IusZZ); |
||
| 1510 | FCF(10,8) = C3*IusZY; |
||
| 1511 | FCF(11,6) = C3*IusXZ; |
||
| 1512 | FCF(11,7) = C3*IusYZ; |
||
| 1513 | FCF(11,8) = C3*(IusYY + IusXX); |
||
| 1514 | |||
| 4925 | parduino | 1515 | /*// block44 |
| 4654 | parduino | 1516 | FCF(9,9) = C1*HstuXX + C3*(HstuYY + HstuZZ); |
| 4925 | parduino | 1517 | FCF(9,10) = C4*HstuXY; |
| 1518 | FCF(9,11) = C4*HstuXZ; |
||
| 1519 | FCF(10,9) = C4*HstuXY; |
||
| 4654 | parduino | 1520 | FCF(10,10) = C1*HstuYY + C3*(HstuXX + HstuZZ); |
| 4925 | parduino | 1521 | FCF(10,11) = C4*HstuYZ; |
| 1522 | FCF(11,9) = C4*HstuXZ; |
||
| 1523 | FCF(11,10) = C4*HstuYZ; |
||
| 1524 | FCF(11,11) = C1*HstuZZ + C3*(HstuYY + HstuXX);*/ |
||
| 4654 | parduino | 1525 | |
| 4925 | parduino | 1526 | // block44 |
| 1527 | FCF(9,9) = C3*(HstuYY + HstuZZ); |
||
| 1528 | FCF(9,10) = C3*HstuXY; |
||
| 1529 | FCF(9,11) = C3*HstuXZ; |
||
| 1530 | FCF(10,9) = FCF(9,10); |
||
| 1531 | FCF(10,10) = C3*(HstuXX + HstuZZ); |
||
| 1532 | FCF(10,11) = C3*HstuYZ; |
||
| 1533 | FCF(11,9) = FCF(9,11); |
||
| 1534 | FCF(11,10) = FCF(10,11); |
||
| 1535 | FCF(11,11) = C3*(HstuYY + HstuXX); |
||
| 1536 | |||
| 4656 | parduino | 1537 | // enhanced strain portion of the stabilization stiffness matrix |
| 1538 | // define the constitutive coefficients |
||
| 1539 | double CssXX = C1*Jinv(0,0)*Jinv(0,0) + C3*(Jinv(0,1)*Jinv(0,1) + Jinv(0,2)*Jinv(0,2)); |
||
| 4925 | parduino | 1540 | double CssXY = C4*Jinv(0,0)*Jinv(0,1); |
| 1541 | double CssXZ = C4*Jinv(0,0)*Jinv(0,2); |
||
| 4656 | parduino | 1542 | double CssYY = C1*Jinv(0,1)*Jinv(0,1) + C3*(Jinv(0,0)*Jinv(0,0) + Jinv(0,2)*Jinv(0,2)); |
| 4925 | parduino | 1543 | double CssYZ = C4*Jinv(0,1)*Jinv(0,2); |
| 4656 | parduino | 1544 | double CssZZ = C1*Jinv(0,2)*Jinv(0,2) + C3*(Jinv(0,0)*Jinv(0,0) + Jinv(0,1)*Jinv(0,1)); |
| 1545 | |||
| 1546 | double CttXX = C1*Jinv(1,0)*Jinv(1,0) + C3*(Jinv(1,1)*Jinv(1,1) + Jinv(1,2)*Jinv(1,2)); |
||
| 4925 | parduino | 1547 | double CttXY = C4*Jinv(1,0)*Jinv(1,1); |
| 1548 | double CttXZ = C4*Jinv(1,0)*Jinv(1,2); |
||
| 4656 | parduino | 1549 | double CttYY = C1*Jinv(1,1)*Jinv(1,1) + C3*(Jinv(1,0)*Jinv(1,0) + Jinv(1,2)*Jinv(1,2)); |
| 4925 | parduino | 1550 | double CttYZ = C4*Jinv(1,1)*Jinv(1,2); |
| 4656 | parduino | 1551 | double CttZZ = C1*Jinv(1,2)*Jinv(1,2) + C3*(Jinv(1,0)*Jinv(1,0) + Jinv(1,1)*Jinv(1,1)); |
| 1552 | |||
| 1553 | double CuuXX = C1*Jinv(2,0)*Jinv(2,0) + C3*(Jinv(2,1)*Jinv(2,1) + Jinv(2,2)*Jinv(2,2)); |
||
| 4925 | parduino | 1554 | double CuuXY = C4*Jinv(2,0)*Jinv(2,1); |
| 1555 | double CuuXZ = C4*Jinv(2,0)*Jinv(2,2); |
||
| 4656 | parduino | 1556 | double CuuYY = C1*Jinv(2,1)*Jinv(2,1) + C3*(Jinv(2,0)*Jinv(2,0) + Jinv(2,2)*Jinv(2,2)); |
| 4925 | parduino | 1557 | double CuuYZ = C4*Jinv(2,1)*Jinv(2,2); |
| 4656 | parduino | 1558 | double CuuZZ = C1*Jinv(2,2)*Jinv(2,2) + C3*(Jinv(2,0)*Jinv(2,0) + Jinv(2,1)*Jinv(2,1)); |
| 1559 | |||
| 1560 | double CstXX = C1*Jinv(0,0)*Jinv(1,0) + C3*(Jinv(0,1)*Jinv(1,1) + Jinv(0,2)*Jinv(1,2)); |
||
| 1561 | double CstXY = C2*Jinv(0,0)*Jinv(1,1) + C3*Jinv(0,1)*Jinv(1,0); |
||
| 1562 | double CstXZ = C2*Jinv(0,0)*Jinv(1,2) + C3*Jinv(0,2)*Jinv(1,0); |
||
| 1563 | double CstYX = C2*Jinv(0,1)*Jinv(1,0) + C3*Jinv(0,0)*Jinv(1,1); |
||
| 1564 | double CstYY = C1*Jinv(0,1)*Jinv(1,1) + C3*(Jinv(0,0)*Jinv(1,0) + Jinv(0,2)*Jinv(1,2)); |
||
| 1565 | double CstYZ = C2*Jinv(0,1)*Jinv(1,2) + C3*Jinv(0,2)*Jinv(1,1); |
||
| 1566 | double CstZX = C2*Jinv(0,2)*Jinv(1,0) + C3*Jinv(0,0)*Jinv(1,2); |
||
| 1567 | double CstZY = C2*Jinv(0,2)*Jinv(1,1) + C3*Jinv(0,1)*Jinv(1,2); |
||
| 1568 | double CstZZ = C1*Jinv(0,2)*Jinv(1,2) + C3*(Jinv(0,0)*Jinv(1,0) + Jinv(0,1)*Jinv(1,1)); |
||
| 1569 | |||
| 1570 | double CsuXX = C1*Jinv(0,0)*Jinv(2,0) + C3*(Jinv(0,1)*Jinv(2,1) + Jinv(0,2)*Jinv(2,2)); |
||
| 1571 | double CsuXY = C2*Jinv(0,0)*Jinv(2,1) + C3*Jinv(0,1)*Jinv(2,0); |
||
| 1572 | double CsuXZ = C2*Jinv(0,0)*Jinv(2,2) + C3*Jinv(0,2)*Jinv(2,0); |
||
| 1573 | double CsuYX = C2*Jinv(0,1)*Jinv(2,0) + C3*Jinv(0,0)*Jinv(2,1); |
||
| 1574 | double CsuYY = C1*Jinv(0,1)*Jinv(2,1) + C3*(Jinv(0,0)*Jinv(2,0) + Jinv(0,2)*Jinv(2,2)); |
||
| 1575 | double CsuYZ = C2*Jinv(0,1)*Jinv(2,2) + C3*Jinv(0,2)*Jinv(2,1); |
||
| 1576 | double CsuZX = C2*Jinv(0,2)*Jinv(2,0) + C3*Jinv(0,0)*Jinv(2,2); |
||
| 1577 | double CsuZY = C2*Jinv(0,2)*Jinv(2,1) + C3*Jinv(0,1)*Jinv(2,2); |
||
| 1578 | double CsuZZ = C1*Jinv(0,2)*Jinv(2,2) + C3*(Jinv(0,0)*Jinv(2,0) + Jinv(0,1)*Jinv(2,1)); |
||
| 1579 | |||
| 1580 | double CtuXX = C1*Jinv(1,0)*Jinv(2,0) + C3*(Jinv(1,1)*Jinv(2,1) + Jinv(1,2)*Jinv(2,2)); |
||
| 1581 | double CtuXY = C2*Jinv(1,0)*Jinv(2,1) + C3*Jinv(1,1)*Jinv(2,0); |
||
| 1582 | double CtuXZ = C2*Jinv(1,0)*Jinv(2,2) + C3*Jinv(1,2)*Jinv(2,0); |
||
| 1583 | double CtuYX = C2*Jinv(1,1)*Jinv(2,0) + C3*Jinv(1,0)*Jinv(2,1); |
||
| 1584 | double CtuYY = C1*Jinv(1,1)*Jinv(2,1) + C3*(Jinv(1,0)*Jinv(2,0) + Jinv(1,2)*Jinv(2,2)); |
||
| 1585 | double CtuYZ = C2*Jinv(1,1)*Jinv(2,2) + C3*Jinv(1,2)*Jinv(2,1); |
||
| 1586 | double CtuZX = C2*Jinv(1,2)*Jinv(2,0) + C3*Jinv(1,0)*Jinv(2,2); |
||
| 1587 | double CtuZY = C2*Jinv(1,2)*Jinv(2,1) + C3*Jinv(1,1)*Jinv(2,2); |
||
| 1588 | double CtuZZ = C1*Jinv(1,2)*Jinv(2,2) + C3*(Jinv(1,0)*Jinv(2,0) + Jinv(1,1)*Jinv(2,1)); |
||
| 1589 | |||
| 1590 | // define the integrated matrix [FenT][C][Fen] |
||
| 1591 | Matrix FeCFe(9,9); |
||
| 1592 | Matrix FeCFeInv(9,9); |
||
| 1593 | |||
| 1594 | // block11 |
||
| 1595 | FeCFe(0,0) = CssXX*J0789; |
||
| 1596 | FeCFe(0,1) = CssXY*J0789; |
||
| 1597 | FeCFe(0,2) = CssXZ*J0789; |
||
| 4925 | parduino | 1598 | FeCFe(1,0) = FeCFe(0,1); |
| 4656 | parduino | 1599 | FeCFe(1,1) = CssYY*J0789; |
| 1600 | FeCFe(1,2) = CssYZ*J0789; |
||
| 4925 | parduino | 1601 | FeCFe(2,0) = FeCFe(0,2); |
| 1602 | FeCFe(2,1) = FeCFe(1,2); |
||
| 4656 | parduino | 1603 | FeCFe(2,2) = CssZZ*J0789; |
| 1604 | |||
| 1605 | // block12 |
||
| 1606 | FeCFe(0,3) = CstXX*J619; |
||
| 1607 | FeCFe(0,4) = CstXY*J619; |
||
| 1608 | FeCFe(0,5) = CstXZ*J619; |
||
| 1609 | FeCFe(1,3) = CstYX*J619; |
||
| 1610 | FeCFe(1,4) = CstYY*J619; |
||
| 1611 | FeCFe(1,5) = CstYZ*J619; |
||
| 1612 | FeCFe(2,3) = CstZX*J619; |
||
| 1613 | FeCFe(2,4) = CstZY*J619; |
||
| 1614 | FeCFe(2,5) = CstZZ*J619; |
||
| 1615 | |||
| 1616 | // block13 |
||
| 1617 | FeCFe(0,6) = CsuXX*J518; |
||
| 1618 | FeCFe(0,7) = CsuXY*J518; |
||
| 1619 | FeCFe(0,8) = CsuXZ*J518; |
||
| 1620 | FeCFe(1,6) = CsuYX*J518; |
||
| 1621 | FeCFe(1,7) = CsuYY*J518; |
||
| 1622 | FeCFe(1,8) = CsuYZ*J518; |
||
| 1623 | FeCFe(2,6) = CsuZX*J518; |
||
| 1624 | FeCFe(2,7) = CsuZY*J518; |
||
| 1625 | FeCFe(2,8) = CsuZZ*J518; |
||
| 1626 | |||
| 1627 | // block21 |
||
| 1628 | FeCFe(3,0) = CstXX*J619; |
||
| 1629 | FeCFe(3,1) = CstYX*J619; |
||
| 1630 | FeCFe(3,2) = CstZX*J619; |
||
| 4925 | parduino | 1631 | FeCFe(4,0) = CstXY*J619; |
| 1632 | FeCFe(4,1) = CstYY*J619; |
||
| 4656 | parduino | 1633 | FeCFe(4,2) = CstZY*J619; |
| 4925 | parduino | 1634 | FeCFe(5,0) = CstXZ*J619; |
| 1635 | FeCFe(5,1) = CstYZ*J619; |
||
| 4656 | parduino | 1636 | FeCFe(5,2) = CstZZ*J619; |
| 1637 | |||
| 1638 | // block22 |
||
| 1639 | FeCFe(3,3) = CttXX*J0879; |
||
| 1640 | FeCFe(3,4) = CttXY*J0879; |
||
| 1641 | FeCFe(3,5) = CttXZ*J0879; |
||
| 4925 | parduino | 1642 | FeCFe(4,3) = FeCFe(3,4); |
| 4656 | parduino | 1643 | FeCFe(4,4) = CttYY*J0879; |
| 1644 | FeCFe(4,5) = CttYZ*J0879; |
||
| 4925 | parduino | 1645 | FeCFe(5,3) = FeCFe(3,5); |
| 1646 | FeCFe(5,4) = FeCFe(4,5); |
||
| 4656 | parduino | 1647 | FeCFe(5,5) = CttZZ*J0879; |
| 1648 | |||
| 1649 | // block23 |
||
| 1650 | FeCFe(3,6) = CtuXX*J417; |
||
| 1651 | FeCFe(3,7) = CtuXY*J417; |
||
| 1652 | FeCFe(3,8) = CtuXZ*J417; |
||
| 1653 | FeCFe(4,6) = CtuYX*J417; |
||
| 1654 | FeCFe(4,7) = CtuYY*J417; |
||
| 1655 | FeCFe(4,8) = CtuYZ*J417; |
||
| 1656 | FeCFe(5,6) = CtuZX*J417; |
||
| 1657 | FeCFe(5,7) = CtuZY*J417; |
||
| 1658 | FeCFe(5,8) = CtuZZ*J417; |
||
| 1659 | |||
| 1660 | // block31 |
||
| 1661 | FeCFe(6,0) = CsuXX*J518; |
||
| 1662 | FeCFe(6,1) = CsuYX*J518; |
||
| 1663 | FeCFe(6,2) = CsuZX*J518; |
||
| 4925 | parduino | 1664 | FeCFe(7,0) = CsuXY*J518; |
| 1665 | FeCFe(7,1) = CsuYY*J518; |
||
| 4656 | parduino | 1666 | FeCFe(7,2) = CsuZY*J518; |
| 4925 | parduino | 1667 | FeCFe(8,0) = CsuXZ*J518; |
| 1668 | FeCFe(8,1) = CsuYZ*J518; |
||
| 4656 | parduino | 1669 | FeCFe(8,2) = CsuZZ*J518; |
| 1670 | |||
| 1671 | // block32 |
||
| 1672 | FeCFe(6,3) = CtuXX*J417; |
||
| 1673 | FeCFe(6,4) = CtuYX*J417; |
||
| 1674 | FeCFe(6,5) = CtuZX*J417; |
||
| 4925 | parduino | 1675 | FeCFe(7,3) = CtuXY*J417; |
| 1676 | FeCFe(7,4) = CtuYY*J417; |
||
| 4656 | parduino | 1677 | FeCFe(7,5) = CtuZY*J417; |
| 4925 | parduino | 1678 | FeCFe(8,3) = CtuXZ*J417; |
| 1679 | FeCFe(8,4) = CtuYZ*J417; |
||
| 4656 | parduino | 1680 | FeCFe(8,5) = CtuZZ*J417; |
| 1681 | |||
| 1682 | // block33 |
||
| 1683 | FeCFe(6,6) = CuuXX*J0978; |
||
| 1684 | FeCFe(6,7) = CuuXY*J0978; |
||
| 1685 | FeCFe(6,8) = CuuXZ*J0978; |
||
| 4925 | parduino | 1686 | FeCFe(7,6) = FeCFe(6,7); |
| 4656 | parduino | 1687 | FeCFe(7,7) = CuuYY*J0978; |
| 1688 | FeCFe(7,8) = CuuYZ*J0978; |
||
| 4925 | parduino | 1689 | FeCFe(8,6) = FeCFe(6,8); |
| 1690 | FeCFe(8,7) = FeCFe(7,8); |
||
| 4656 | parduino | 1691 | FeCFe(8,8) = CuuZZ*J0978; |
| 1692 | |||
| 1693 | // inverse of [FenT][C][Fen] |
||
| 1694 | FeCFe.Invert(FeCFeInv); |
||
| 1695 | |||
| 1696 | // define the integrated matrix [FenT][C][Fhg] |
||
| 1697 | Matrix FeCFhg(9,12); |
||
| 1698 | FeCFhg.Zero(); |
||
| 1699 | |||
| 1700 | // block11 |
||
| 1701 | FeCFhg(0,0) = CstXX*J0789 + CssXX*J619; |
||
| 1702 | FeCFhg(0,1) = CstXY*J0789 + CssXY*J619; |
||
| 1703 | FeCFhg(0,2) = CstXZ*J0789 + CssXZ*J619; |
||
| 1704 | FeCFhg(1,0) = CstYX*J0789 + CssXY*J619; |
||
| 1705 | FeCFhg(1,1) = CstYY*J0789 + CssYY*J619; |
||
| 1706 | FeCFhg(1,2) = CstYZ*J0789 + CssYZ*J619; |
||
| 1707 | FeCFhg(2,0) = CstZX*J0789 + CssXZ*J619; |
||
| 1708 | FeCFhg(2,1) = CstZY*J0789 + CssYZ*J619; |
||
| 1709 | FeCFhg(2,2) = CstZZ*J0789 + CssZZ*J619; |
||
| 1710 | |||
| 1711 | // block12 |
||
| 1712 | FeCFhg(0,3) = CsuXX*J619 + CstXX*J518; |
||
| 1713 | FeCFhg(0,4) = CsuXY*J619 + CstXY*J518; |
||
| 1714 | FeCFhg(0,5) = CsuXZ*J619 + CstXZ*J518; |
||
| 1715 | FeCFhg(1,3) = CsuYX*J619 + CstYX*J518; |
||
| 1716 | FeCFhg(1,4) = CsuYY*J619 + CstYY*J518; |
||
| 1717 | FeCFhg(1,5) = CsuYZ*J619 + CstYZ*J518; |
||
| 1718 | FeCFhg(2,3) = CsuZX*J619 + CstZX*J518; |
||
| 1719 | FeCFhg(2,4) = CsuZY*J619 + CstZY*J518; |
||
| 1720 | FeCFhg(2,5) = CsuZZ*J619 + CstZZ*J518; |
||
| 1721 | |||
| 1722 | // block13 |
||
| 1723 | FeCFhg(0,6) = CsuXX*J0789 + CssXX*J518; |
||
| 1724 | FeCFhg(0,7) = CsuXY*J0789 + CssXY*J518; |
||
| 1725 | FeCFhg(0,8) = CsuXZ*J0789 + CssXZ*J518; |
||
| 1726 | FeCFhg(1,6) = CsuYX*J0789 + CssXY*J518; |
||
| 1727 | FeCFhg(1,7) = CsuYY*J0789 + CssYY*J518; |
||
| 1728 | FeCFhg(1,8) = CsuYZ*J0789 + CssYZ*J518; |
||
| 1729 | FeCFhg(2,6) = CsuZX*J0789 + CssXZ*J518; |
||
| 1730 | FeCFhg(2,7) = CsuZY*J0789 + CssYZ*J518; |
||
| 1731 | FeCFhg(2,8) = CsuZZ*J0789 + CssZZ*J518; |
||
| 1732 | |||
| 1733 | // block21 |
||
| 1734 | FeCFhg(3,0) = CstXX*J0879 + CttXX*J619; |
||
| 1735 | FeCFhg(3,1) = CstYX*J0879 + CttXY*J619; |
||
| 1736 | FeCFhg(3,2) = CstZX*J0879 + CttXZ*J619; |
||
| 1737 | FeCFhg(4,0) = CstXY*J0879 + CttXY*J619; |
||
| 1738 | FeCFhg(4,1) = CstYY*J0879 + CttYY*J619; |
||
| 1739 | FeCFhg(4,2) = CstZY*J0879 + CttYZ*J619; |
||
| 1740 | FeCFhg(5,0) = CstXZ*J0879 + CttXZ*J619; |
||
| 1741 | FeCFhg(5,1) = CstYZ*J0879 + CttYZ*J619; |
||
| 1742 | FeCFhg(5,2) = CstZZ*J0879 + CttZZ*J619; |
||
| 1743 | |||
| 1744 | // block22 |
||
| 1745 | FeCFhg(3,3) = CtuXX*J0879 + CttXX*J417; |
||
| 1746 | FeCFhg(3,4) = CtuXY*J0879 + CttXY*J417; |
||
| 1747 | FeCFhg(3,5) = CtuXZ*J0879 + CttXZ*J417; |
||
| 1748 | FeCFhg(4,3) = CtuYX*J0879 + CttXY*J417; |
||
| 1749 | FeCFhg(4,4) = CtuYY*J0879 + CttYY*J417; |
||
| 1750 | FeCFhg(4,5) = CtuYZ*J0879 + CttYZ*J417; |
||
| 1751 | FeCFhg(5,3) = CtuZX*J0879 + CttXZ*J417; |
||
| 1752 | FeCFhg(5,4) = CtuZY*J0879 + CttYZ*J417; |
||
| 1753 | FeCFhg(5,5) = CtuZZ*J0879 + CttZZ*J417; |
||
| 1754 | |||
| 1755 | // block23 |
||
| 1756 | FeCFhg(3,6) = CtuXX*J619 + CstXX*J417; |
||
| 1757 | FeCFhg(3,7) = CtuXY*J619 + CstYX*J417; |
||
| 1758 | FeCFhg(3,8) = CtuXZ*J619 + CstZX*J417; |
||
| 1759 | FeCFhg(4,6) = CtuYX*J619 + CstXY*J417; |
||
| 1760 | FeCFhg(4,7) = CtuYY*J619 + CstYY*J417; |
||
| 1761 | FeCFhg(4,8) = CtuYZ*J619 + CstZY*J417; |
||
| 1762 | FeCFhg(5,6) = CtuZX*J619 + CstXZ*J417; |
||
| 1763 | FeCFhg(5,7) = CtuZY*J619 + CstYZ*J417; |
||
| 1764 | FeCFhg(5,8) = CtuZZ*J619 + CstZZ*J417; |
||
| 1765 | |||
| 1766 | // block31 |
||
| 1767 | FeCFhg(6,0) = CtuXX*J518 + CsuXX*J417; |
||
| 1768 | FeCFhg(6,1) = CtuYX*J518 + CsuYX*J417; |
||
| 1769 | FeCFhg(6,2) = CtuZX*J518 + CsuZX*J417; |
||
| 1770 | FeCFhg(7,0) = CtuXY*J518 + CsuXY*J417; |
||
| 1771 | FeCFhg(7,1) = CtuYY*J518 + CsuYY*J417; |
||
| 1772 | FeCFhg(7,2) = CtuZY*J518 + CsuZY*J417; |
||
| 1773 | FeCFhg(8,0) = CtuXZ*J518 + CsuXZ*J417; |
||
| 1774 | FeCFhg(8,1) = CtuYZ*J518 + CsuYZ*J417; |
||
| 1775 | FeCFhg(8,2) = CtuZZ*J518 + CsuZZ*J417; |
||
| 1776 | |||
| 1777 | // block32 |
||
| 1778 | FeCFhg(6,3) = CtuXX*J0978 + CuuXX*J417; |
||
| 1779 | FeCFhg(6,4) = CtuYX*J0978 + CuuXY*J417; |
||
| 1780 | FeCFhg(6,5) = CtuZX*J0978 + CuuXZ*J417; |
||
| 1781 | FeCFhg(7,3) = CtuXY*J0978 + CuuXY*J417; |
||
| 1782 | FeCFhg(7,4) = CtuYY*J0978 + CuuYY*J417; |
||
| 1783 | FeCFhg(7,5) = CtuZY*J0978 + CuuYZ*J417; |
||
| 1784 | FeCFhg(8,3) = CtuXZ*J0978 + CuuXZ*J417; |
||
| 1785 | FeCFhg(8,4) = CtuYZ*J0978 + CuuYZ*J417; |
||
| 1786 | FeCFhg(8,5) = CtuZZ*J0978 + CuuZZ*J417; |
||
| 1787 | |||
| 1788 | // block33 |
||
| 1789 | FeCFhg(6,6) = CsuXX*J0978 + CuuXX*J518; |
||
| 1790 | FeCFhg(6,7) = CsuYX*J0978 + CuuXY*J518; |
||
| 1791 | FeCFhg(6,8) = CsuZX*J0978 + CuuXZ*J518; |
||
| 1792 | FeCFhg(7,6) = CsuXY*J0978 + CuuXY*J518; |
||
| 1793 | FeCFhg(7,7) = CsuYY*J0978 + CuuYY*J518; |
||
| 1794 | FeCFhg(7,8) = CsuZY*J0978 + CuuYZ*J518; |
||
| 1795 | FeCFhg(8,6) = CsuXZ*J0978 + CuuXZ*J518; |
||
| 1796 | FeCFhg(8,7) = CsuYZ*J0978 + CuuYZ*J518; |
||
| 1797 | FeCFhg(8,8) = CsuZZ*J0978 + CuuZZ*J518; |
||
| 1798 | |||
| 4925 | parduino | 1799 | /*// off-diagonal enhanced strain matrix [Fe]^T[C][F] is same as [Fe]^T[C][Fhg] with exception of last three columns |
| 1800 | Matrix FeCF(9,12); |
||
| 1801 | FeCF = FeCFhg; |
||
| 4656 | parduino | 1802 | |
| 4925 | parduino | 1803 | // block 14 |
| 1804 | FeCF(0,9) = C3*(J21411*(Jinv(2,1)*Jinv(0,1) + Jinv(2,2)*Jinv(0,2)) + J31310*(Jinv(1,1)*Jinv(0,1) + Jinv(1,2)*Jinv(0,2)) + J16*(Jinv(0,1)*Jinv(0,1) + Jinv(0,2)*Jinv(0,2))); |
||
| 1805 | FeCF(0,10) = C3*(J21411*Jinv(2,0)*Jinv(0,1) + J31310*Jinv(1,0)*Jinv(0,1) + J16*Jinv(0,0)*Jinv(0,1)); |
||
| 1806 | FeCF(0,11) = C3*(J21411*Jinv(2,0)*Jinv(0,2) + J31310*Jinv(1,0)*Jinv(0,2) + J16*Jinv(0,0)*Jinv(0,2)); |
||
| 1807 | FeCF(1,9) = C3*(J21411*Jinv(2,1)*Jinv(0,0) + J31310*Jinv(1,1)*Jinv(0,0) + J16*Jinv(0,0)*Jinv(0,1)); |
||
| 1808 | FeCF(1,10) = C3*(J21411*(Jinv(2,0)*Jinv(0,0) + Jinv(2,2)*Jinv(0,2)) + J31310*(Jinv(1,0)*Jinv(0,0) + Jinv(1,2)*Jinv(0,2)) + J16*(Jinv(0,0)*Jinv(0,0) + Jinv(0,2)*Jinv(0,2))); |
||
| 1809 | FeCF(1,11) = C3*(J21411*Jinv(2,1)*Jinv(0,2) + J31310*Jinv(1,1)*Jinv(0,2) + J16*Jinv(0,1)*Jinv(0,2)); |
||
| 1810 | FeCF(2,9) = C3*(J21411*Jinv(2,2)*Jinv(0,0) + J31310*Jinv(1,2)*Jinv(0,0) + J16*Jinv(0,0)*Jinv(0,2)); |
||
| 1811 | FeCF(2,10) = C3*(J21411*Jinv(2,2)*Jinv(0,1) + J31310*Jinv(1,2)*Jinv(0,1) + J16*Jinv(0,1)*Jinv(0,2)); |
||
| 1812 | FeCF(2,11) = C3*(J21411*(Jinv(2,0)*Jinv(0,0) + Jinv(2,1)*Jinv(0,1)) + J31310*(Jinv(1,0)*Jinv(0,0) + Jinv(1,1)*Jinv(0,1)) + J16*(Jinv(0,0)*Jinv(0,0) + Jinv(0,1)*Jinv(0,1))); |
||
| 4656 | parduino | 1813 | |
| 4925 | parduino | 1814 | // block 24 |
| 1815 | FeCF(3,9) = C3*(J11512*(Jinv(2,1)*Jinv(1,1) + Jinv(2,2)*Jinv(1,2)) + J31013*(Jinv(1,1)*Jinv(0,1) + Jinv(1,2)*Jinv(0,2)) + J16*(Jinv(1,1)*Jinv(1,1) + Jinv(1,2)*Jinv(1,2))); |
||
| 1816 | FeCF(3,10) = C3*(J11512*Jinv(2,0)*Jinv(1,1) + J31013*Jinv(1,1)*Jinv(0,0) + J16*Jinv(1,0)*Jinv(1,1)); |
||
| 1817 | FeCF(3,11) = C3*(J11512*Jinv(2,0)*Jinv(1,2) + J31013*Jinv(1,2)*Jinv(0,0) + J16*Jinv(1,0)*Jinv(1,2)); |
||
| 1818 | FeCF(4,9) = C3*(J11512*Jinv(2,1)*Jinv(1,0) + J31013*Jinv(1,0)*Jinv(0,1) + J16*Jinv(1,0)*Jinv(1,1)); |
||
| 1819 | FeCF(4,10) = C3*(J11512*(Jinv(2,0)*Jinv(1,0) + Jinv(2,2)*Jinv(1,2)) + J31013*(Jinv(1,0)*Jinv(0,0) + Jinv(1,2)*Jinv(0,2)) + J16*(Jinv(1,0)*Jinv(1,0) + Jinv(1,2)*Jinv(1,2))); |
||
| 1820 | FeCF(4,11) = C3*(J11512*Jinv(2,1)*Jinv(1,2) + J31013*Jinv(1,2)*Jinv(0,1) + J16*Jinv(1,1)*Jinv(1,2)); |
||
| 1821 | FeCF(5,9) = C3*(J11512*Jinv(2,2)*Jinv(1,0) + J31013*Jinv(1,0)*Jinv(0,2) + J16*Jinv(1,0)*Jinv(1,2)); |
||
| 1822 | FeCF(5,10) = C3*(J11512*Jinv(2,2)*Jinv(1,1) + J31013*Jinv(1,1)*Jinv(0,2) + J16*Jinv(1,1)*Jinv(1,2)); |
||
| 1823 | FeCF(5,11) = C3*(J11512*(Jinv(2,0)*Jinv(1,0) + Jinv(2,1)*Jinv(1,1)) + J31013*(Jinv(1,0)*Jinv(0,0) + Jinv(1,1)*Jinv(0,1)) + J16*(Jinv(1,0)*Jinv(1,0) + Jinv(1,1)*Jinv(1,1))); |
||
| 1824 | |||
| 1825 | // block 34 |
||
| 1826 | FeCF(6,9) = C3*(J11215*(Jinv(2,1)*Jinv(1,1) + Jinv(2,2)*Jinv(1,2)) + J21114*(Jinv(2,1)*Jinv(0,1) + Jinv(2,2)*Jinv(0,2)) + J16*(Jinv(2,1)*Jinv(2,1) + Jinv(2,2)*Jinv(2,2))); |
||
| 1827 | FeCF(6,10) = C3*(J11215*Jinv(2,1)*Jinv(1,0) + J21114*Jinv(2,1)*Jinv(0,0) + J16*Jinv(2,0)*Jinv(2,1)); |
||
| 1828 | FeCF(6,11) = C3*(J11215*Jinv(2,2)*Jinv(1,0) + J21114*Jinv(2,2)*Jinv(0,0) + J16*Jinv(2,0)*Jinv(2,2)); |
||
| 1829 | FeCF(7,9) = C3*(J11215*Jinv(2,0)*Jinv(1,1) + J21114*Jinv(2,0)*Jinv(0,1) + J16*Jinv(2,0)*Jinv(2,1)); |
||
| 1830 | FeCF(7,10) = C3*(J11215*(Jinv(2,0)*Jinv(1,0) + Jinv(2,2)*Jinv(1,2)) + J21114*(Jinv(2,0)*Jinv(0,0) + Jinv(2,2)*Jinv(0,2)) + J16*(Jinv(2,0)*Jinv(2,0) + Jinv(2,2)*Jinv(2,2))); |
||
| 1831 | FeCF(7,11) = C3*(J11215*Jinv(2,2)*Jinv(1,1) + J21114*Jinv(2,2)*Jinv(0,1) + J16*Jinv(2,1)*Jinv(2,2)); |
||
| 1832 | FeCF(8,9) = C3*(J11215*Jinv(2,0)*Jinv(1,2) + J21114*Jinv(2,0)*Jinv(0,2) + J16*Jinv(2,0)*Jinv(2,2)); |
||
| 1833 | FeCF(8,10) = C3*(J11215*Jinv(2,1)*Jinv(1,2) + J21114*Jinv(2,1)*Jinv(0,2) + J16*Jinv(2,1)*Jinv(2,2)); |
||
| 1834 | FeCF(8,11) = C3*(J11215*(Jinv(2,0)*Jinv(1,0) + Jinv(2,1)*Jinv(1,1)) + J21114*(Jinv(2,0)*Jinv(0,0) + Jinv(2,1)*Jinv(0,1)) + J16*(Jinv(2,0)*Jinv(2,0) + Jinv(2,1)*Jinv(2,1))); |
||
| 1835 | |||
| 1836 | // transpose the FeCF matrix to get the FCFe matrix |
||
| 1837 | Matrix FCFe(12,9); |
||
| 1838 | FCFe = Transpose(9,12,FeCF);*/ |
||
| 1839 | |||
| 1840 | // transpose the Kwu matrix |
||
| 1841 | Matrix KuT(12,9); |
||
| 1842 | KuT = Transpose(9,12,FeCFhg); |
||
| 1843 | |||
| 1844 | // compute the stabilization stiffness matrix |
||
| 1845 | Kstab.Zero(); |
||
| 1846 | |||
| 1847 | /*Matrix KKF(12,12); |
||
| 1848 | Matrix FKK(12,12); |
||
| 1849 | Matrix KKFKK(12,12); |
||
| 1850 | |||
| 1851 | KKF = KuT*FeCFeInv*FeCF; |
||
| 1852 | |