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