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