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:307­319 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