1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/shelleas.h,v 1.13 2017/01/12 14:46:44 masarati Exp $ */ 2 /* 3 * MBDyn (C) is a multibody analysis code. 4 * http://www.mbdyn.org 5 * 6 * Copyright (C) 2010-2017 7 * 8 * Marco Morandini <morandini@aero.polimi.it> 9 * Riccardo Vescovini <vescovini@aero.polimi.it> 10 * 11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 12 * via La Masa, 34 - 20156 Milano, Italy 13 * http://www.aero.polimi.it 14 * 15 * Changing this copyright notice is forbidden. 16 * 17 * This program is free software; you can redistribute it and/or modify 18 * it under the terms of the GNU General Public License as published by 19 * the Free Software Foundation (version 2 of the License). 20 * 21 * 22 * This program is distributed in the hope that it will be useful, 23 * but WITHOUT ANY WARRANTY; without even the implied warranty of 24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 * GNU General Public License for more details. 26 * 27 * You should have received a copy of the GNU General Public License 28 * along with this program; if not, write to the Free Software 29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 30 */ 31 32 /* 33 * Inspired by 34 * Wojciech Witkowski 35 * "4-Node combined shell element with semi-EAS-ANS strain interpolations 36 * in 6-parameter shell theories with drilling degrees of freedom" 37 * Comput Mech (2009) 43:307319 DOI 10.1007/s00466-008-0307-x 38 */ 39 40 #ifndef SHELLEAS_H 41 #define SHELLEAS_H 42 43 #include "myassert.h" 44 #include "except.h" 45 46 #include "strnode.h" 47 #include "elem.h" 48 49 #include "constltp.h" 50 51 /* da spostare in .cc */ 52 #include "Rot.hh" 53 #include "joint.h" 54 55 #include "shell.h" 56 57 // Shell4EAS - begin 58 59 class Shell4EAS 60 : virtual public Elem, 61 public Shell 62 { 63 #if 0 64 protected: 65 static const unsigned int iNumPrivData = 66 +3 // 0 ( 1 -> 3) - strain 67 +3 // 3 ( 4 -> 6) - curvature 68 +3 // 6 ( 7 -> 9) - force 69 +3 // 9 (10 -> 12) - moment 70 +3 // 12 (13 -> 15) - position 71 +3 // 15 (16 -> 18) - orientation vector 72 +3 // 18 (19 -> 21) - angular velocity 73 +3 // 21 (22 -> 24) - strain rate 74 +3 // 24 (25 -> 27) - curvature rate 75 ; 76 77 static unsigned int iGetPrivDataIdx_int(const char *s, 78 ConstLawType::Type type); 79 #endif 80 private: 81 public: 82 // numbered according to 83 // 84 // ^ 85 // 4 o-----+-----o 3 86 // | 1_2 | 2_2 | 87 // --+-----+-----+-> 88 // | 1_1 | 2_1 | 89 // 1 o-----+-----o 2 90 // 91 enum IntegrationPoint { 92 IP_1_1 = 0, 93 IP_1_2 = 1, 94 IP_1_3 = 2, 95 IP_2_1 = 3, 96 IP_2_2 = 4, 97 IP_2_3 = 5, 98 IP_3_1 = 6, 99 IP_3_2 = 7, 100 IP_3_3 = 8, 101 102 NUMIP = 4 103 }; 104 105 static doublereal xi_i[NUMIP][2]; 106 static doublereal w_i[NUMIP]; 107 108 // numbered according to the side they are defined on 109 enum ShearStrainEvaluationPoint { 110 SSEP_1 = 0, 111 SSEP_2 = 1, 112 SSEP_3 = 2, 113 SSEP_4 = 3, 114 }; 115 116 enum NodeName { 117 NODE1 = 0, 118 NODE2 = 1, 119 NODE3 = 2, 120 NODE4 = 3, 121 122 NUMNODES = 4 123 }; 124 125 static doublereal xi_n[NUMNODES][2]; 126 127 static doublereal xi_0[2]; 128 129 enum Deformations { 130 STRAIN = 0, 131 CURVAT = 1, 132 133 NUMDEFORM = 2 134 }; 135 136 protected: 137 // Pointers to nodes 138 const StructNode* pNode[NUMNODES]; 139 140 #if 0 141 // Node offsets - TODO: offsets 142 const Vec3 f[NUMNODES]; 143 Vec3 fRef[NUMNODES]; 144 #endif 145 146 // nodal positions (0: initial; otherwise current) 147 Vec3 xa_0[NUMNODES]; 148 Vec3 xa[NUMNODES]; 149 // current nodal orientation 150 Mat3x3 iTa[NUMNODES]; 151 Mat3x3 iTa_i[NUMIP]; 152 // Euler vector of Ra 153 Vec3 phi_tilde_n[NUMNODES]; 154 155 // Average orientation matrix 156 Vec3 phi_tilde_i[NUMIP]; 157 // Average orientation matrix 158 // .. in reference configuration 159 Mat3x3 T0_overline; 160 // .. in current configuration 161 Mat3x3 T_overline; 162 163 // Orientation matrix 164 // .. in reference configuration 165 Mat3x3 T_0_0; 166 Mat3x3 T_0_i[NUMIP]; 167 // .. in current configuration 168 Mat3x3 T_0; 169 Mat3x3 T_i[NUMIP]; 170 171 Mat3x3 Phi_Delta_i[NUMIP][NUMNODES]; 172 Mat3x3 Kappa_delta_i_1[NUMIP][NUMNODES]; 173 Mat3x3 Kappa_delta_i_2[NUMIP][NUMNODES]; 174 175 // // Orientation matrix of the shear strain evaluation points 176 // Mat3x3 R[NUMSSEP]; 177 // Mat3x3 RRef[NUMSSEP]; 178 // Mat3x3 RPrev[NUMSSEP]; 179 180 // rotation tensors 181 Mat3x3 Q_i[NUMIP]; 182 183 // Orientation tensor derivative axial vector 184 Vec3 k_1_i[NUMIP]; 185 Vec3 k_2_i[NUMIP]; 186 // Mat3x3 T_1_i[NUMIP]; 187 // Mat3x3 T_2_i[NUMIP]; 188 189 // linear deformation vectors 190 // .. in reference configuration 191 Vec3 eps_tilde_1_0_i[NUMIP]; 192 Vec3 eps_tilde_2_0_i[NUMIP]; 193 // .. in current configuration 194 Vec3 eps_tilde_1_i[NUMIP]; 195 Vec3 eps_tilde_2_i[NUMIP]; 196 197 // angular deformation vectors 198 // .. in reference configuration 199 Vec3 k_tilde_1_0_i[NUMIP]; 200 Vec3 k_tilde_2_0_i[NUMIP]; 201 // .. in current configuration 202 Vec3 k_tilde_1_i[NUMIP]; 203 Vec3 k_tilde_2_i[NUMIP]; 204 205 206 #if 0 207 // Angular velocity of the sections - TODO: viscoelastic 208 Vec3 Omega[NUMSEZ]; 209 Vec3 OmegaRef[NUMSEZ]; 210 #endif 211 212 // Temporary data - TODO 213 #if 0 214 Vec6 Az[NUMSEZ]; 215 Vec6 AzRef[NUMSEZ]; 216 Vec6 AzLoc[NUMSEZ]; 217 Vec6 DefLoc[NUMSEZ]; 218 Vec6 DefLocRef[NUMSEZ]; 219 Vec6 DefLocPrev[NUMSEZ]; 220 #endif 221 222 protected: 223 fmh S_alpha_beta_0; 224 doublereal alpha_0; 225 doublereal alpha_i[NUMIP]; 226 vfmh L_alpha_beta_i; 227 228 vfmh B_overline_i; 229 230 vfmh P_i; 231 232 Vec3 y_i_1[NUMIP]; 233 Vec3 y_i_2[NUMIP]; 234 235 vh beta; 236 vh epsilon_hat; 237 vh epsilon; 238 239 #ifdef USE_CL_IN_SHELL 240 // Constitutive law handlers at integration points 241 ConstitutiveLawOwner<vh, fmh>* pD[NUMIP]; 242 #else // ! USE_CL_IN_SHELL 243 vvh FRef; 244 bool bPreStress; 245 vh PreStress; 246 #endif // ! USE_CL_IN_SHELL 247 248 // Reference constitutive law tangent matrices 249 vfmh DRef; 250 251 //stress 252 vvh stress_i; 253 254 // Is first residual 255 bool bFirstRes; 256 257 // for derived elements that add external contributions 258 // to internal forces 259 virtual void AddInternalForces(Vec6 &,unsigned int)260 AddInternalForces(Vec6& /* AzLoc */ , unsigned int /* iSez */ ) { 261 NO_OP; 262 }; 263 264 private: 265 void UpdateNodalAndAveragePosAndOrientation(); 266 void ComputeInitialNodeAndIptOrientation(); 267 void InterpolateOrientation(); 268 // void ComputeIPSEPRotations(); 269 void ComputeIPCurvature(); 270 271 public: 272 Shell4EAS(unsigned int uL, 273 const DofOwner* pDO, 274 const StructNode* pN1, const StructNode* pN2, 275 const StructNode* pN3, const StructNode* pN4, 276 #if 0 // TODO: offset 277 const Vec3& f1, const Vec3& f2, 278 const Vec3& f3, const Vec3& f4, 279 #endif 280 const Mat3x3& R1, const Mat3x3& R2, 281 const Mat3x3& R3, const Mat3x3& R4, 282 #ifdef USE_CL_IN_SHELL 283 const ConstitutiveLaw<vh, fmh>** pDTmp, 284 #else // ! USE_CL_IN_SHELL 285 const fmh& pDTmp, 286 const vh& PreStress, 287 #endif // ! USE_CL_IN_SHELL 288 flag fOut); 289 290 virtual ~Shell4EAS(void); 291 292 // Shell type GetShellType(void)293 virtual Shell::Type GetShellType(void) const { 294 return Shell::ELASTIC; 295 }; 296 iGetNumDof(void)297 virtual unsigned int iGetNumDof(void) const { 298 //return 14; 299 return 13; 300 }; 301 GetDofType(unsigned int i)302 virtual DofOrder::Order GetDofType(unsigned int i) const { 303 ASSERT(i >= 0 && i < iGetNumDof()); 304 return DofOrder::ALGEBRAIC; 305 }; 306 GetEqType(unsigned int i)307 virtual DofOrder::Order GetEqType(unsigned int i) const { 308 ASSERT(i >= 0 && i < iGetNumDof()); 309 return DofOrder::ALGEBRAIC; 310 }; 311 312 // Contribution to restart file 313 virtual std::ostream& Restart(std::ostream& out) const; 314 315 #if 0 316 virtual void 317 AfterConvergence(const VectorHandler& X, const VectorHandler& XP); 318 319 // Inverse dynamics 320 virtual void 321 AfterConvergence(const VectorHandler& X, const VectorHandler& XP, 322 const VectorHandler& XPP); 323 #endif 324 325 // Workspace size 326 virtual void WorkSpaceDim(integer * piNumRows,integer * piNumCols)327 WorkSpaceDim(integer* piNumRows, integer* piNumCols) const { 328 *piNumRows = 6*4 + iGetNumDof(); 329 *piNumCols = 6*4 + iGetNumDof(); 330 }; 331 332 // Initial settings 333 void 334 SetValue(DataManager *pDM, 335 VectorHandler& /* X */ , VectorHandler& /* XP */ , 336 SimulationEntity::Hints *ph = 0); 337 338 #if 0 339 // Prepares reference parameters after prediction 340 virtual void 341 AfterPredict(VectorHandler& /* X */ , VectorHandler& /* XP */ ); 342 #endif 343 344 // Residual assembly 345 virtual SubVectorHandler& AssRes(SubVectorHandler& WorkVec, 346 doublereal dCoef, 347 const VectorHandler& XCurr, 348 const VectorHandler& XPrimeCurr); 349 350 #if 0 351 // Inverse dynamics 352 virtual SubVectorHandler& 353 AssRes(SubVectorHandler& WorkVec, 354 const VectorHandler& XCurr , 355 const VectorHandler& XPrimeCurr , 356 const VectorHandler& XPrimePrimeCurr , 357 InverseDynamics::Order iOrder = InverseDynamics::INVERSE_DYNAMICS); 358 #endif 359 360 // Jacobian matrix assembly 361 virtual VariableSubMatrixHandler& 362 AssJac(VariableSubMatrixHandler& WorkMat, 363 doublereal dCoef, 364 const VectorHandler& XCurr, 365 const VectorHandler& XPrimeCurr); 366 367 #if 0 368 // Matrix assembly for eigenvalues 369 void 370 AssMats(VariableSubMatrixHandler& WorkMatA, 371 VariableSubMatrixHandler& WorkMatB, 372 const VectorHandler& XCurr, 373 const VectorHandler& XPrimeCurr); 374 #endif 375 376 virtual void Output(OutputHandler& OH) const; 377 iGetInitialNumDof(void)378 virtual unsigned int iGetInitialNumDof(void) const { 379 return 0; 380 }; 381 382 virtual void InitialWorkSpaceDim(integer * piNumRows,integer * piNumCols)383 InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const { 384 *piNumRows = 6*4; 385 *piNumCols = 6*4; 386 }; 387 SetInitialValue(VectorHandler &)388 virtual void SetInitialValue(VectorHandler& /* X */ ) { 389 NO_OP; 390 }; 391 392 virtual VariableSubMatrixHandler& 393 InitialAssJac(VariableSubMatrixHandler& WorkMat, 394 const VectorHandler& XCurr); 395 396 virtual SubVectorHandler& 397 InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr); 398 399 #if 0 400 // Access to private data 401 virtual unsigned int iGetNumPrivData(void) const; 402 virtual unsigned int iGetPrivDataIdx(const char *s) const; 403 virtual doublereal dGetPrivData(unsigned int i) const; 404 #endif 405 406 // Access to nodes 407 virtual const StructNode* pGetNode(unsigned int i) const; 408 409 /******** PER IL SOLUTORE PARALLELO *********/ 410 /* Fornisce il tipo e la label dei nodi che sono connessi all'elemento 411 * utile per l'assemblaggio della matrice di connessione fra i dofs */ 412 virtual void GetConnectedNodes(std::vector<const Node * > & connectedNodes)413 GetConnectedNodes(std::vector<const Node *>& connectedNodes) const { 414 connectedNodes.resize(NUMNODES); 415 for (int i = 0; i < NUMNODES; i++) { 416 connectedNodes[i] = pNode[i]; 417 } 418 }; 419 /**************************************************/ 420 }; 421 422 // Shell4EAS - end 423 424 #endif // SHELLEAS_H 425 426