1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora, Radu Serban
13 // =============================================================================
14 // Four node shell with geom.exact kinematics
15 // =============================================================================
16
17 #include "chrono/core/ChException.h"
18 #include "chrono/physics/ChSystem.h"
19 #include "chrono/timestepper/ChState.h"
20 #include "chrono/fea/ChElementShellReissner4.h"
21 #include "chrono/fea/ChRotUtils.h"
22 #include <cmath>
23
24 #define CHUSE_ANS
25 //#define CHUSE_EAS
26 #define CHUSE_KGEOMETRIC
27 //#define CHSIMPLIFY_DROT
28
29 namespace chrono {
30 namespace fea {
31
32 //--------------------------------------------------------------
33 // utility functions
34 //--------------------------------------------------------------
35
L1(const double xi[2])36 static inline double L1(const double xi[2]) {
37 return 0.25 * (1. + xi[0]) * (1. + xi[1]);
38 };
39
L2(const double xi[2])40 static inline double L2(const double xi[2]) {
41 return 0.25 * (1. - xi[0]) * (1. + xi[1]);
42 };
43
L3(const double xi[2])44 static inline double L3(const double xi[2]) {
45 return 0.25 * (1. - xi[0]) * (1. - xi[1]);
46 };
47
L4(const double xi[2])48 static inline double L4(const double xi[2]) {
49 return 0.25 * (1. + xi[0]) * (1. - xi[1]);
50 };
51
52 typedef double (*LI_Type)(const double xi[2]);
53 static LI_Type LI[4] = {&L1, &L2, &L3, &L4};
54
L1_1(const double xi[2])55 static inline double L1_1(const double xi[2]) {
56 return 0.25 * (1. + xi[1]);
57 }
58
L1_2(const double xi[2])59 static inline double L1_2(const double xi[2]) {
60 return 0.25 * (1. + xi[0]);
61 }
62
L2_1(const double xi[2])63 static inline double L2_1(const double xi[2]) {
64 return -0.25 * (1. + xi[1]);
65 }
66
L2_2(const double xi[2])67 static inline double L2_2(const double xi[2]) {
68 return 0.25 * (1. - xi[0]);
69 }
70
L3_1(const double xi[2])71 static inline double L3_1(const double xi[2]) {
72 return -0.25 * (1. - xi[1]);
73 }
74
L3_2(const double xi[2])75 static inline double L3_2(const double xi[2]) {
76 return -0.25 * (1. - xi[0]);
77 }
78
L4_1(const double xi[2])79 static inline double L4_1(const double xi[2]) {
80 return 0.25 * (1. - xi[1]);
81 }
82
L4_2(const double xi[2])83 static inline double L4_2(const double xi[2]) {
84 return -0.25 * (1. + xi[0]);
85 }
86
87 typedef double (*LI_J_Type)(const double xi[2]);
88 static LI_J_Type LI_J[4][2] = {
89 {&L1_1, &L1_2},
90 {&L2_1, &L2_2},
91 {&L3_1, &L3_2},
92 {&L4_1, &L4_2},
93 };
94
Interp(const ChVector<> * const v,const double xi[2])95 static ChVector<> Interp(const ChVector<>* const v, const double xi[2]) {
96 ChVector<> r = v[0] * L1(xi) + v[1] * L2(xi) + v[2] * L3(xi) + v[3] * L4(xi);
97 return r;
98 }
99
InterpDeriv1(const ChVector<> * const v,const ChMatrixNM<double,4,2> & der_mat)100 static ChVector<> InterpDeriv1(const ChVector<>* const v, const ChMatrixNM<double, 4, 2>& der_mat) {
101 ChVector<> r = v[0] * der_mat(0, 0) + v[1] * der_mat(1, 0) + v[2] * der_mat(2, 0) + v[3] * der_mat(3, 0);
102 return r;
103 }
104
InterpDeriv2(const ChVector<> * const v,const ChMatrixNM<double,4,2> & der_mat)105 static ChVector<> InterpDeriv2(const ChVector<>* const v, const ChMatrixNM<double, 4, 2>& der_mat) {
106 ChVector<> r = v[0] * der_mat(0, 1) + v[1] * der_mat(1, 1) + v[2] * der_mat(2, 1) + v[3] * der_mat(3, 1);
107 return r;
108 }
109
InterpDeriv(const ChVector<> * const v,const ChMatrixNM<double,4,2> & der_mat,ChVector<> & der1,ChVector<> & der2)110 static void InterpDeriv(const ChVector<>* const v,
111 const ChMatrixNM<double, 4, 2>& der_mat,
112 ChVector<>& der1,
113 ChVector<>& der2) {
114 der1 = InterpDeriv1(v, der_mat);
115 der2 = InterpDeriv2(v, der_mat);
116 return;
117 }
118
InterpDeriv_xi1(const ChVector<> * const v,const double xi[2])119 static ChVector<> InterpDeriv_xi1(const ChVector<>* const v, const double xi[2]) {
120 ChVector<> r = v[0] * L1_1(xi) + v[1] * L2_1(xi) + v[2] * L3_1(xi) + v[3] * L4_1(xi);
121 return r;
122 }
123
InterpDeriv_xi2(const ChVector<> * const v,const double xi[2])124 static ChVector<> InterpDeriv_xi2(const ChVector<>* const v, const double xi[2]) {
125 ChVector<> r = v[0] * L1_2(xi) + v[1] * L2_2(xi) + v[2] * L3_2(xi) + v[3] * L4_2(xi);
126 return r;
127 }
128
129 //--------------------------------------------------------------
130 // ChElementShellReissner4
131 //--------------------------------------------------------------
132
133 // Static integration tables
134
135 double ChElementShellReissner4::xi_i[ChElementShellReissner4::NUMIP][2] = {{-1. / std::sqrt(3.), -1. / std::sqrt(3.)},
136 {1. / std::sqrt(3.), -1. / std::sqrt(3.)},
137 {1. / std::sqrt(3.), 1. / std::sqrt(3.)},
138 {-1. / std::sqrt(3.), 1. / std::sqrt(3.)}};
139
140 double ChElementShellReissner4::w_i[ChElementShellReissner4::NUMIP] = {1., 1., 1., 1.};
141
142 double ChElementShellReissner4::xi_A[ChElementShellReissner4::NUMSSEP][2] = {{0., 1.}, {-1., 0.}, {0., -1.}, {1., 0.}};
143
144 double ChElementShellReissner4::xi_n[ChElementShellReissner4::NUMNODES][2] = {{1., 1.},
145 {-1., 1.},
146 {-1., -1.},
147 {1., -1.}};
148
149 double ChElementShellReissner4::xi_0[2] = {0., 0.};
150
UpdateNodalAndAveragePosAndOrientation()151 void ChElementShellReissner4::UpdateNodalAndAveragePosAndOrientation() {
152 ChMatrix33<> Tn[NUMNODES];
153 ChMatrix33<> T_avg;
154 T_avg.setZero();
155 for (int i = 0; i < NUMNODES; i++) {
156 xa[i] = this->m_nodes[i]->GetPos();
157 Tn[i] = this->m_nodes[i]->GetA() * iTa[i];
158 T_avg += this->m_nodes[i]->GetA() * iTa[i]; //***TODO*** use predicted rot?
159 }
160 T_avg *= 0.25;
161 T_overline = rotutils::Rot(rotutils::VecRot(T_avg));
162 /*
163 // ***ALEX*** test an alternative for T_overline:
164 // average four rotations with quaternion averaging:
165 ChQuaternion<> qTa = Tn[0].Get_A_quaternion();
166 ChQuaternion<> qTb = Tn[1].Get_A_quaternion();
167 if ( (qTa ^ qTb) < 0)
168 qTb *= -1;
169 ChQuaternion<> qTc = Tn[2].Get_A_quaternion();
170 if ( (qTa ^ qTc) < 0)
171 qTc *= -1;
172 ChQuaternion<> qTd = Tn[3].Get_A_quaternion();
173 if ( (qTa ^ qTd) < 0)
174 qTd *= -1;
175 ChQuaternion<> Tavg =(qTa + qTb + qTc + qTd ).GetNormalized();
176 T_overline.Set_A_quaternion(Tavg);
177 */
178
179 ChMatrix33<> R_tilde_n[NUMNODES];
180 for (int i = 0; i < NUMNODES; i++) {
181 R_tilde_n[i] = T_overline.transpose() * Tn[i];
182 phi_tilde_n[i] = rotutils::VecRot(R_tilde_n[i]);
183 // if (phi_tilde_n[i].Length()*CH_C_RAD_TO_DEG > 15)
184 // GetLog() << "WARNING phi_tilde_n[" << i << "]=" << phi_tilde_n[i].Length()*CH_C_RAD_TO_DEG << "�\n";
185 }
186 }
187
ComputeInitialNodeOrientation()188 void ChElementShellReissner4::ComputeInitialNodeOrientation() {
189 for (int i = 0; i < NUMNODES; i++) {
190 xa[i] = this->m_nodes[i]->GetPos();
191 }
192 for (int i = 0; i < NUMNODES; i++) {
193 ChVector<> t1 = InterpDeriv_xi1(xa, xi_n[i]);
194 t1 = t1 / t1.Length();
195 ChVector<> t2 = InterpDeriv_xi2(xa, xi_n[i]);
196 t2 = t2 / t2.Length();
197 ChVector<> t3 = Vcross(t1, t2);
198 t3 = t3 / t3.Length();
199 t2 = Vcross(t3, t1);
200 ChMatrix33<> t123(t1, t2, t3);
201 iTa[i] = m_nodes[i]->GetA().transpose() * t123;
202 }
203 for (int i = 0; i < NUMIP; i++) {
204 iTa_i[i] = ChMatrix33<>(1);
205 }
206 for (int i = 0; i < NUMSSEP; i++) {
207 iTa_A[i] = ChMatrix33<>(1);
208 }
209 UpdateNodalAndAveragePosAndOrientation();
210 InterpolateOrientation();
211 for (int i = 0; i < NUMIP; i++) {
212 ChVector<> t1 = InterpDeriv_xi1(xa, xi_i[i]);
213 t1 = t1 / t1.Length();
214 ChVector<> t2 = InterpDeriv_xi2(xa, xi_i[i]);
215 t2 = t2 / t2.Length();
216 ChVector<> t3 = Vcross(t1, t2);
217 t3 = t3 / t3.Length();
218 t2 = Vcross(t3, t1);
219 ChMatrix33<> t123(t1, t2, t3);
220 iTa_i[i] = T_i[i].transpose() * t123;
221 }
222 for (int i = 0; i < NUMSSEP; i++) {
223 ChVector<> t1 = InterpDeriv_xi1(xa, xi_A[i]);
224 t1 = t1 / t1.Length();
225 ChVector<> t2 = InterpDeriv_xi2(xa, xi_A[i]);
226 t2 = t2 / t2.Length();
227 ChVector<> t3 = Vcross(t1, t2);
228 t3 = t3 / t3.Length();
229 t2 = Vcross(t3, t1);
230 ChMatrix33<> t123(t1, t2, t3);
231 iTa_A[i] = T_A[i].transpose() * t123;
232 }
233 InterpolateOrientation();
234 }
235
InterpolateOrientation()236 void ChElementShellReissner4::InterpolateOrientation() {
237 ChMatrix33<> DRot_I_phi_tilde_n_MT_T_overline[NUMNODES];
238 ChMatrix33<> Ri, Gammai;
239 for (int n = 0; n < NUMNODES; n++) {
240 ChMatrix33<> mDrot_I = rotutils::DRot_I(phi_tilde_n[n]);
241 #ifdef CHSIMPLIFY_DROT
242 mDrot_I = ChMatrix33<>(1);
243 #endif
244 DRot_I_phi_tilde_n_MT_T_overline[n] = mDrot_I * T_overline.transpose();
245 }
246 for (int i = 0; i < NUMIP; i++) {
247 phi_tilde_i[i] = Interp(phi_tilde_n, xi_i[i]);
248 rotutils::RotAndDRot(phi_tilde_i[i], Ri, Gammai);
249 #ifdef CHSIMPLIFY_DROT
250 Gammai = ChMatrix33<>(1);
251 #endif
252 T_i[i] = T_overline * Ri * iTa_i[i];
253 ChMatrix33<> T_overline_Gamma_tilde_i(T_overline * Gammai);
254 for (int n = 0; n < NUMNODES; n++) {
255 Phi_Delta_i[i][n] = T_overline_Gamma_tilde_i * DRot_I_phi_tilde_n_MT_T_overline[n];
256 }
257 }
258 ChVector<> phi_tilde_0 = Interp(phi_tilde_n, xi_0);
259 T_0 = T_overline * rotutils::Rot(phi_tilde_0);
260 for (int i = 0; i < NUMSSEP; i++) {
261 phi_tilde_A[i] = Interp(phi_tilde_n, xi_A[i]);
262 rotutils::RotAndDRot(phi_tilde_A[i], Ri, Gammai);
263 #ifdef CHSIMPLIFY_DROT
264 Gammai = ChMatrix33<>(1);
265 #endif
266 T_A[i] = T_overline * Ri * iTa_A[i];
267 ChMatrix33<> T_overline_Gamma_tilde_A(T_overline * Gammai);
268 for (int n = 0; n < NUMNODES; n++) {
269 Phi_Delta_A[i][n] = T_overline_Gamma_tilde_A * DRot_I_phi_tilde_n_MT_T_overline[n];
270 }
271 }
272 }
273
ComputeIPCurvature()274 void ChElementShellReissner4::ComputeIPCurvature() {
275 ChMatrix33<> Gamma_I_n_MT_T_overline[NUMNODES];
276 for (int n = 0; n < NUMNODES; n++) {
277 ChMatrix33<> mDrot_I = rotutils::DRot_I(phi_tilde_n[n]);
278 #ifdef CHSIMPLIFY_DROT
279 mDrot_I = ChMatrix33<>(1);
280 #endif
281 Gamma_I_n_MT_T_overline[n] = mDrot_I * T_overline.transpose();
282 }
283 for (int i = 0; i < NUMIP; i++) {
284 ChVector<> phi_tilde_1_i;
285 ChVector<> phi_tilde_2_i;
286 InterpDeriv(phi_tilde_n, L_alpha_beta_i[i], phi_tilde_1_i, phi_tilde_2_i);
287 ChMatrix33<> mGamma_tilde_i = rotutils::DRot(phi_tilde_i[i]);
288 #ifdef CHSIMPLIFY_DROT
289 mGamma_tilde_i = ChMatrix33<>(1);
290 #endif
291 ChMatrix33<> T_overlineGamma_tilde_i(T_overline * mGamma_tilde_i);
292 k_1_i[i] = T_overlineGamma_tilde_i * phi_tilde_1_i;
293 k_2_i[i] = T_overlineGamma_tilde_i * phi_tilde_2_i;
294 ChMatrix33<> tmp1 = T_overline * rotutils::Elle(phi_tilde_i[i], phi_tilde_1_i);
295 ChMatrix33<> tmp2 = T_overline * rotutils::Elle(phi_tilde_i[i], phi_tilde_2_i);
296 #ifdef CHSIMPLIFY_DROT
297 tmp1 = T_overline;
298 tmp2 = T_overline;
299 #endif
300 for (int n = 0; n < NUMNODES; n++) {
301 Kappa_delta_i_1[i][n] =
302 tmp1 * Gamma_I_n_MT_T_overline[n] * LI[n](xi_i[i]) + Phi_Delta_i[i][n] * L_alpha_beta_i[i](n, 0);
303 Kappa_delta_i_2[i][n] =
304 tmp2 * Gamma_I_n_MT_T_overline[n] * LI[n](xi_i[i]) + Phi_Delta_i[i][n] * L_alpha_beta_i[i](n, 1);
305 }
306 }
307 }
308
ChElementShellReissner4()309 ChElementShellReissner4::ChElementShellReissner4() : tot_thickness(0) {
310 m_nodes.resize(4);
311
312 // add here other constructor-time initializations
313 for (int i = 0; i < NUMIP; i++) {
314 L_alpha_beta_i[i].setZero();
315 B_overline_i[i].setZero();
316 D_overline_i[i].setZero();
317 G_i[i].setZero();
318 P_i[i].setZero();
319 }
320
321 for (int i = 0; i < NUMSSEP; i++)
322 L_alpha_beta_A[i].setZero();
323 }
324
~ChElementShellReissner4()325 ChElementShellReissner4::~ChElementShellReissner4() {}
326
SetNodes(std::shared_ptr<ChNodeFEAxyzrot> nodeA,std::shared_ptr<ChNodeFEAxyzrot> nodeB,std::shared_ptr<ChNodeFEAxyzrot> nodeC,std::shared_ptr<ChNodeFEAxyzrot> nodeD)327 void ChElementShellReissner4::SetNodes(std::shared_ptr<ChNodeFEAxyzrot> nodeA,
328 std::shared_ptr<ChNodeFEAxyzrot> nodeB,
329 std::shared_ptr<ChNodeFEAxyzrot> nodeC,
330 std::shared_ptr<ChNodeFEAxyzrot> nodeD) {
331 assert(nodeA);
332 assert(nodeB);
333 assert(nodeC);
334 assert(nodeD);
335
336 m_nodes[0] = nodeA;
337 m_nodes[1] = nodeB;
338 m_nodes[2] = nodeC;
339 m_nodes[3] = nodeD;
340 std::vector<ChVariables*> mvars;
341 mvars.push_back(&m_nodes[0]->Variables());
342 mvars.push_back(&m_nodes[1]->Variables());
343 mvars.push_back(&m_nodes[2]->Variables());
344 mvars.push_back(&m_nodes[3]->Variables());
345 Kmatr.SetVariables(mvars);
346 }
347
EvaluateGP(int igp)348 ChVector<> ChElementShellReissner4::EvaluateGP(int igp) {
349 return GetNodeA()->GetPos() * L1(xi_i[igp]) + GetNodeB()->GetPos() * L2(xi_i[igp]) +
350 GetNodeC()->GetPos() * L3(xi_i[igp]) + GetNodeD()->GetPos() * L4(xi_i[igp]);
351 }
EvaluatePT(int ipt)352 ChVector<> ChElementShellReissner4::EvaluatePT(int ipt) {
353 return GetNodeA()->GetPos() * L1(xi_n[ipt]) + GetNodeB()->GetPos() * L2(xi_n[ipt]) +
354 GetNodeC()->GetPos() * L3(xi_n[ipt]) + GetNodeD()->GetPos() * L4(xi_n[ipt]);
355 }
356
357 // -----------------------------------------------------------------------------
358 // Add a layer.
359 // -----------------------------------------------------------------------------
360
AddLayer(double thickness,double theta,std::shared_ptr<ChMaterialShellReissner> material)361 void ChElementShellReissner4::AddLayer(double thickness,
362 double theta,
363 std::shared_ptr<ChMaterialShellReissner> material) {
364 m_layers.push_back(Layer(this, thickness, theta, material));
365 SetLayerZreferenceCentered();
366 }
367
SetLayerZreferenceCentered()368 void ChElementShellReissner4::SetLayerZreferenceCentered() {
369 // accumulate element thickness.
370 tot_thickness = 0;
371 for (size_t kl = 0; kl < m_layers.size(); kl++) {
372 tot_thickness += m_layers[kl].Get_thickness();
373 }
374
375 // Loop again over the layers and calculate the z levels of layers, by centering them
376 m_layers_z.clear();
377 m_layers_z.push_back(-0.5 * this->GetThickness());
378 for (size_t kl = 0; kl < m_layers.size(); kl++) {
379 m_layers_z.push_back(m_layers_z[kl] + m_layers[kl].Get_thickness());
380 }
381 }
382
SetLayerZreference(double z_from_bottom)383 void ChElementShellReissner4::SetLayerZreference(double z_from_bottom) {
384 // accumulate element thickness.
385 tot_thickness = 0;
386 for (size_t kl = 0; kl < m_layers.size(); kl++) {
387 tot_thickness += m_layers[kl].Get_thickness();
388 }
389
390 // Loop again over the layers and calculate the z levels of layers, by centering them
391 m_layers_z.clear();
392 m_layers_z.push_back(z_from_bottom);
393 for (size_t kl = 0; kl < m_layers.size(); kl++) {
394 m_layers_z.push_back(m_layers_z[kl] + m_layers[kl].Get_thickness());
395 }
396 }
397
398 // -----------------------------------------------------------------------------
399 // Set as neutral position
400 // -----------------------------------------------------------------------------
401
402 // Initial element setup.
SetAsNeutral()403 void ChElementShellReissner4::SetAsNeutral() {
404 GetNodeA()->GetX0ref().SetPos(GetNodeA()->GetPos());
405 GetNodeB()->GetX0ref().SetPos(GetNodeB()->GetPos());
406 GetNodeC()->GetX0ref().SetPos(GetNodeC()->GetPos());
407 GetNodeD()->GetX0ref().SetPos(GetNodeD()->GetPos());
408 GetNodeA()->GetX0ref().SetRot(GetNodeA()->GetRot());
409 GetNodeB()->GetX0ref().SetRot(GetNodeB()->GetRot());
410 GetNodeC()->GetX0ref().SetRot(GetNodeC()->GetRot());
411 GetNodeD()->GetX0ref().SetRot(GetNodeD()->GetRot());
412 }
413
414 // -----------------------------------------------------------------------------
415 // Interface to ChElementBase base class
416 // -----------------------------------------------------------------------------
417
418 // Initial element setup.
SetupInitial(ChSystem * system)419 void ChElementShellReissner4::SetupInitial(ChSystem* system) {
420 // Align initial pos/rot of nodes to actual pos/rot
421 SetAsNeutral();
422
423 ComputeInitialNodeOrientation();
424 // UpdateNodalAndAveragePosAndOrientation();
425 // InterpolateOrientation();
426 // copy ref values
427 T0_overline = T_overline;
428 T_0_0 = T_0;
429 for (int i = 0; i < NUMNODES; i++) {
430 xa_0[i] = xa[i];
431 }
432 for (int i = 0; i < NUMIP; i++) {
433 T_0_i[i] = T_i[i];
434 }
435 for (int i = 0; i < NUMSSEP; i++) {
436 T_0_A[i] = T_A[i];
437 }
438
439 ChMatrixNM<double, 4, 4> M_0;
440 ChMatrixNM<double, 4, 4> M_0_Inv;
441 {
442 ChVector<> x_1 = InterpDeriv_xi1(xa, xi_0);
443 ChVector<> x_2 = InterpDeriv_xi2(xa, xi_0);
444 S_alpha_beta_0(0, 0) = T_0_0.Get_A_Xaxis() ^ x_1;
445 S_alpha_beta_0(1, 0) = T_0_0.Get_A_Yaxis() ^ x_1;
446 S_alpha_beta_0(0, 1) = T_0_0.Get_A_Xaxis() ^ x_2;
447 S_alpha_beta_0(1, 1) = T_0_0.Get_A_Yaxis() ^ x_2;
448 alpha_0 = S_alpha_beta_0(0, 0) * S_alpha_beta_0(1, 1) - S_alpha_beta_0(0, 1) * S_alpha_beta_0(1, 0);
449
450 M_0(0, 0) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(0, 0);
451 M_0(0, 1) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(0, 1);
452 M_0(0, 2) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(0, 0);
453 M_0(0, 3) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(0, 1);
454
455 M_0(1, 0) = S_alpha_beta_0(1, 0) * S_alpha_beta_0(1, 0);
456 M_0(1, 1) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(1, 1);
457 M_0(1, 2) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(1, 0);
458 M_0(1, 3) = S_alpha_beta_0(1, 0) * S_alpha_beta_0(1, 1);
459
460 M_0(2, 0) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(1, 0);
461 M_0(2, 1) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(1, 1);
462 M_0(2, 2) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(1, 1);
463 M_0(2, 3) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(1, 0);
464
465 M_0(3, 0) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(1, 0);
466 M_0(3, 1) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(1, 1);
467 M_0(3, 2) = S_alpha_beta_0(0, 1) * S_alpha_beta_0(1, 0);
468 M_0(3, 3) = S_alpha_beta_0(0, 0) * S_alpha_beta_0(1, 1);
469
470 M_0_Inv = M_0.inverse();
471 }
472
473 for (int i = 0; i < NUMIP; i++) {
474 ChMatrixNM<double, 4, 2> L_alpha_B_i;
475 ChVector<> x_1 = InterpDeriv_xi1(xa, xi_i[i]);
476 ChVector<> x_2 = InterpDeriv_xi2(xa, xi_i[i]);
477 S_alpha_beta_i[i](0, 0) = T_0_i[i].Get_A_Xaxis() ^ x_1;
478 S_alpha_beta_i[i](1, 0) = T_0_i[i].Get_A_Yaxis() ^ x_1;
479 S_alpha_beta_i[i](0, 1) = T_0_i[i].Get_A_Xaxis() ^ x_2;
480 S_alpha_beta_i[i](1, 1) = T_0_i[i].Get_A_Yaxis() ^ x_2;
481 // alpha_i = det(S_alpha_beta_i)
482 alpha_i[i] =
483 S_alpha_beta_i[i](0, 0) * S_alpha_beta_i[i](1, 1) - S_alpha_beta_i[i](0, 1) * S_alpha_beta_i[i](1, 0);
484
485 ChMatrixNM<double, 2, 2> xi_i_i;
486 xi_i_i = S_alpha_beta_i[i].inverse();
487
488 for (int n = 0; n < NUMNODES; n++) {
489 for (int ii = 0; ii < 2; ii++) {
490 L_alpha_B_i(n, ii) = LI_J[n][ii](xi_i[i]);
491 }
492 }
493
494 L_alpha_beta_i[i] = L_alpha_B_i * xi_i_i;
495
496 double t = xi_i[i][0] * xi_i[i][1];
497 ChMatrixNM<double, 4, IDOFS> H;
498 H.setZero();
499
500 H(0, 0) = xi_i[i][0];
501 H(0, 1) = t;
502
503 H(1, 2) = xi_i[i][1];
504 H(1, 3) = t;
505
506 H(2, 4) = xi_i[i][0];
507 H(2, 5) = t;
508
509 H(3, 6) = xi_i[i][1];
510 H(3, 5) = t;
511
512 ChMatrixNM<double, 12, 4> Perm;
513 Perm.setZero();
514
515 // 1, 5, 4, 2, 3, 6
516 Perm(0, 0) = 1.;
517 Perm(1, 2) = 1.;
518 Perm(3, 3) = 1.;
519 Perm(4, 1) = 1.;
520
521 ChMatrixNM<double, 4, IDOFS> tmpP = M_0_Inv.transpose() * H;
522
523 P_i[i] = (alpha_0 / alpha_i[i]) * Perm * tmpP;
524 }
525 // save initial axial values
526 ComputeIPCurvature();
527 for (int i = 0; i < NUMIP; i++) {
528 InterpDeriv(xa, L_alpha_beta_i[i], y_i_1[i], y_i_2[i]);
529 eps_tilde_1_0_i[i] = T_i[i].transpose() * y_i_1[i];
530 eps_tilde_2_0_i[i] = T_i[i].transpose() * y_i_2[i];
531 k_tilde_1_0_i[i] = T_i[i].transpose() * k_1_i[i];
532 k_tilde_2_0_i[i] = T_i[i].transpose() * k_2_i[i];
533 }
534 for (int i = 0; i < NUMSSEP; i++) {
535 ChMatrixNM<double, 4, 2> L_alpha_B_A;
536 ChVector<> x_1 = InterpDeriv_xi1(xa, xi_A[i]);
537 ChVector<> x_2 = InterpDeriv_xi2(xa, xi_A[i]);
538 S_alpha_beta_A[i](0, 0) = T_0_A[i].Get_A_Xaxis() ^ x_1;
539 S_alpha_beta_A[i](1, 0) = T_0_A[i].Get_A_Yaxis() ^ x_1;
540 S_alpha_beta_A[i](0, 1) = T_0_A[i].Get_A_Xaxis() ^ x_2;
541 S_alpha_beta_A[i](1, 1) = T_0_A[i].Get_A_Yaxis() ^ x_2;
542
543 ChVector<> y_A_1;
544 ChVector<> y_A_2;
545 InterpDeriv(xa, L_alpha_beta_A[i], y_A_1, y_A_2);
546 eps_tilde_1_0_A[i] = T_A[i].transpose() * y_A_1;
547 eps_tilde_2_0_A[i] = T_A[i].transpose() * y_A_2;
548
549 // xi_A_i = S_alpha_beta_A^{-1}
550 ChMatrixNM<double, 2, 2> xi_A_i;
551 xi_A_i = S_alpha_beta_A[i].inverse();
552
553 for (int n = 0; n < NUMNODES; n++) {
554 for (int ii = 0; ii < 2; ii++) {
555 L_alpha_B_A(n, ii) = LI_J[n][ii](xi_A[i]);
556 }
557 }
558
559 L_alpha_beta_A[i] = L_alpha_B_A * xi_A_i;
560 }
561 {
562 ChVector<> y_A_1;
563 ChVector<> y_A_2;
564 for (int i = 0; i < NUMSSEP; i++) {
565 InterpDeriv(xa, L_alpha_beta_A[i], y_A_1, y_A_2);
566 eps_tilde_1_0_A[i] = T_A[i].transpose() * y_A_1;
567 eps_tilde_2_0_A[i] = T_A[i].transpose() * y_A_2;
568 }
569 }
570
571 // Perform layer initialization
572 for (size_t kl = 0; kl < m_layers.size(); kl++) {
573 m_layers[kl].SetupInitial();
574 }
575
576 // compute initial sizes (just for auxiliary information)
577 m_lenX =
578 (0.5 * (GetNodeA()->coord.pos + GetNodeD()->coord.pos) - 0.5 * (GetNodeB()->coord.pos + GetNodeC()->coord.pos))
579 .Length();
580 m_lenY =
581 (0.5 * (GetNodeA()->coord.pos + GetNodeB()->coord.pos) - 0.5 * (GetNodeD()->coord.pos + GetNodeC()->coord.pos))
582 .Length();
583
584 // Compute mass matrix
585 ComputeMassMatrix();
586 }
587
588 // State update.
Update()589 void ChElementShellReissner4::Update() {
590 ChElementGeneric::Update();
591 }
592
593 // Fill the D vector with the current field values at the element nodes.
GetStateBlock(ChVectorDynamic<> & mD)594 void ChElementShellReissner4::GetStateBlock(ChVectorDynamic<>& mD) {
595 mD.resize(4 * 7, 1);
596
597 mD.segment(0, 3) = m_nodes[0]->GetPos().eigen();
598 mD.segment(3, 4) = m_nodes[0]->GetRot().eigen();
599
600 mD.segment(7, 3) = m_nodes[1]->GetPos().eigen();
601 mD.segment(10, 4) = m_nodes[1]->GetRot().eigen();
602
603 mD.segment(14, 3) = m_nodes[2]->GetPos().eigen();
604 mD.segment(17, 4) = m_nodes[2]->GetRot().eigen();
605
606 mD.segment(21, 3) = m_nodes[3]->GetPos().eigen();
607 mD.segment(24, 4) = m_nodes[3]->GetRot().eigen();
608 }
609
610 // Calculate the global matrix H as a linear combination of K, R, and M:
611 // H = Mfactor * [M] + Kfactor * [K] + Rfactor * [R]
612 // NOTE! we assume that this function is computed after one computed
613 // ComputeInternalForces(), that updates inner data for the given node states.
614
ComputeKRMmatricesGlobal(ChMatrixRef H,double Kfactor,double Rfactor,double Mfactor)615 void ChElementShellReissner4::ComputeKRMmatricesGlobal(ChMatrixRef H, double Kfactor, double Rfactor, double Mfactor) {
616 assert((H.rows() == 24) && (H.cols() == 24));
617
618 // Calculate the mass matrix
619 ComputeMassMatrix();
620
621 // Calculate the linear combination Kfactor*[K] + Rfactor*[R]
622 ComputeInternalJacobians(Kfactor, Rfactor);
623
624 // Load Jac + Mfactor*[M] into H
625 for (int i = 0; i < 24; i++)
626 for (int j = 0; j < 24; j++)
627 H(i, j) = m_JacobianMatrix(i, j) + Mfactor * m_MassMatrix(i, j);
628 }
629
630 // Return the mass matrix.
ComputeMmatrixGlobal(ChMatrixRef M)631 void ChElementShellReissner4::ComputeMmatrixGlobal(ChMatrixRef M) {
632 // Calculate the mass matrix
633 ComputeMassMatrix();
634
635 M = m_MassMatrix;
636 }
637
638 // -----------------------------------------------------------------------------
639 // Mass matrix calculation
640 // -----------------------------------------------------------------------------
641
ComputeMassMatrix()642 void ChElementShellReissner4::ComputeMassMatrix() {
643 m_MassMatrix.setZero();
644
645 // loop over all layers, to compute total "mass per area" = sum(rho_i*thickness_i) = average_rho * sum(thickness_i)
646 double avg_density = this->GetDensity();
647 double mass_per_area = avg_density * tot_thickness;
648
649 // Heuristic, simplified 'lumped' mass matrix.
650 // Split the mass in 4 pieces, weighting as the jacobian at integration point, but
651 // lump at the node closest to integration point.
652 // This is simpler than the stiffness-consistent mass matrix that would require
653 // integration over gauss points.
654
655 for (int n = 0; n < NUMNODES; n++) {
656 int igp = (n + 2) % 4; // id of closest gauss point to node
657 double jacobian =
658 alpha_i[igp]; // weight: jacobian at gauss point (another,brute force,option would be 0.25 for all nodes)
659
660 double nodemass = (jacobian * mass_per_area);
661
662 // Approximate (!) inertia of a quarter of tile, note *(1/4) because only the quarter tile,
663 // in local system of Gauss point
664 double Ixx = (pow(this->GetLengthY(), 2) + pow(tot_thickness, 2)) * (1. / 12.) * (1. / 4.) * nodemass;
665 double Iyy = (pow(this->GetLengthX(), 2) + pow(tot_thickness, 2)) * (1. / 12.) * (1. / 4.) * nodemass;
666 double Izz = (pow(this->GetLengthX(), 2) + pow(this->GetLengthY(), 2)) * (1. / 12.) * (1. / 4.) * nodemass;
667 ChMatrix33<> box_inertia(ChVector<>(Ixx, Iyy, Izz));
668 // ..and rotate inertia in local system of node: (local because ystem-level rotational coords of nodes are
669 // ang.vel in loc sys)
670 // I' = A'* T * I * T' * A
671 // = Rot * I * Rot'
672 ChMatrix33<> Rot = m_nodes[n]->GetA().transpose() * T_i[igp];
673 ChMatrix33<> Rot_I(Rot * box_inertia);
674 ChMatrix33<> inertia_n;
675 inertia_n = Rot_I * Rot.transpose();
676
677 int node_off = n * 6;
678 m_MassMatrix(node_off + 0, node_off + 0) = nodemass;
679 m_MassMatrix(node_off + 1, node_off + 1) = nodemass;
680 m_MassMatrix(node_off + 2, node_off + 2) = nodemass;
681 m_MassMatrix.block(node_off + 3, node_off + 3, 3, 3) = inertia_n;
682
683 } // end loop on nodes
684 }
685
686 // -----------------------------------------------------------------------------
687 // Elastic force calculation
688 // -----------------------------------------------------------------------------
689
ComputeInternalForces(ChVectorDynamic<> & Fi)690 void ChElementShellReissner4::ComputeInternalForces(ChVectorDynamic<>& Fi) {
691 Fi.setZero();
692
693 UpdateNodalAndAveragePosAndOrientation();
694 InterpolateOrientation();
695
696 /*
697 for (unsigned int i = 1; i <= iGetNumDof(); i++) {
698 beta(i) = XCurr(iFirstReactionIndex + i);
699 }
700 */ //***TODO*** EAS internal variables not yet implemented
701
702 ComputeIPCurvature();
703 for (int i = 0; i < NUMIP; i++) {
704 InterpDeriv(xa, L_alpha_beta_i[i], y_i_1[i], y_i_2[i]);
705 eps_tilde_1_i[i] = T_i[i].transpose() * y_i_1[i] - eps_tilde_1_0_i[i];
706 eps_tilde_2_i[i] = T_i[i].transpose() * y_i_2[i] - eps_tilde_2_0_i[i];
707 k_tilde_1_i[i] = T_i[i].transpose() * k_1_i[i] - k_tilde_1_0_i[i];
708 k_tilde_2_i[i] = T_i[i].transpose() * k_2_i[i] - k_tilde_2_0_i[i];
709
710 ChStarMatrix33<> myi_1_X(y_i_1[i]);
711 ChStarMatrix33<> myi_2_X(y_i_2[i]);
712 ChStarMatrix33<> mk_1_X(k_1_i[i]);
713 ChStarMatrix33<> mk_2_X(k_2_i[i]);
714 ChMatrix33<> block;
715
716 // parte variabile di B_overline_i
717 for (int n = 0; n < NUMNODES; n++) {
718 ChMatrix33<> Phi_Delta_i_n_LI_i = Phi_Delta_i[i][n] * LI[n](xi_i[i]);
719
720 // delta epsilon_tilde_1_i
721 block = T_i[i] * L_alpha_beta_i[i](n, 0);
722 B_overline_i[i].block(0, 6 * n, 3, 3) = block.transpose();
723 block = T_i[i].transpose() * myi_1_X * Phi_Delta_i_n_LI_i;
724 block = block * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
725 B_overline_i[i].block(0, 3 + 6 * n, 3, 3) = block;
726
727 // delta epsilon_tilde_2_i
728 block = T_i[i] * L_alpha_beta_i[i](n, 1);
729 B_overline_i[i].block(3, 6 * n, 3, 3) = block.transpose();
730 block = T_i[i].transpose() * myi_2_X * Phi_Delta_i_n_LI_i;
731 block = block * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
732 B_overline_i[i].block(3, 3 + 6 * n, 3, 3) = block;
733
734 ChVector<> phi_tilde_1_i;
735 ChVector<> phi_tilde_2_i;
736 InterpDeriv(phi_tilde_n, L_alpha_beta_i[i], phi_tilde_1_i, phi_tilde_2_i);
737
738 // delta k_tilde_1_i
739 block = T_i[i].transpose() * mk_1_X * Phi_Delta_i_n_LI_i + T_i[i].transpose() * Kappa_delta_i_1[i][n];
740 block = block * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
741 B_overline_i[i].block(6, 3 + 6 * n, 3, 3) = block;
742
743 // delta k_tilde_2_i
744 block = T_i[i].transpose() * mk_2_X * Phi_Delta_i_n_LI_i + T_i[i].transpose() * Kappa_delta_i_2[i][n];
745 block = block * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
746 B_overline_i[i].block(9, 3 + 6 * n, 3, 3) = block;
747
748 // delta y_alpha_1
749 block = ChMatrix33<>::Identity() * L_alpha_beta_i[i](n, 0);
750 D_overline_i[i].block(0, 6 * n, 3, 3) = block;
751
752 // delta y_alpha_2
753 block = ChMatrix33<>::Identity() * L_alpha_beta_i[i](n, 1);
754 D_overline_i[i].block(3, 6 * n, 3, 3) = block;
755
756 // delta k_1_i
757 block = Kappa_delta_i_1[i][n] * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
758 D_overline_i[i].block(6, 3 + 6 * n, 3, 3) = block;
759
760 // delta k_2_i
761 block = Kappa_delta_i_2[i][n] * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
762 D_overline_i[i].block(9, 3 + 6 * n, 3, 3) = block;
763
764 // phi_delta
765 block = Phi_Delta_i_n_LI_i * m_nodes[n]->GetA(); //***NEEDED because rotations are body-relative
766 D_overline_i[i].block(12, 3 + 6 * n, 3, 3) = block;
767 }
768 }
769
770 // ANS
771 #ifdef CHUSE_ANS
772
773 ChMatrixNM<double, 6, 24> B_overline_A;
774 ChVector<> y_A_1;
775 ChVector<> y_A_2;
776 ChMatrixNM<double, 4, 24> B_overline_3_ABCD;
777 ChMatrixNM<double, 4, 24> B_overline_6_ABCD;
778
779 for (int i = 0; i < NUMSSEP; i++) {
780 B_overline_A.setZero();
781 InterpDeriv(xa, L_alpha_beta_A[i], y_A_1, y_A_2);
782 eps_tilde_1_A[i] = T_A[i].transpose() * y_A_1 - eps_tilde_1_0_A[i];
783 eps_tilde_2_A[i] = T_A[i].transpose() * y_A_2 - eps_tilde_2_0_A[i];
784
785 ChStarMatrix33<> myA_1_X(y_A_1);
786 ChStarMatrix33<> myA_2_X(y_A_2);
787 ChMatrix33<> block;
788
789 for (int n = 0; n < NUMNODES; n++) {
790 ChMatrix33<> Phi_Delta_A_n_LI_i = Phi_Delta_A[i][n] * LI[n](xi_A[i]);
791
792 // delta epsilon_tilde_1_A
793 block = T_A[i].transpose() * L_alpha_beta_A[i](n, 0);
794 B_overline_A.block(0, 6 * n, 3, 3) = block;
795 block = T_A[i].transpose() * myA_1_X * Phi_Delta_A_n_LI_i;
796 block = block * this->m_nodes[n]->GetA(); //***NEEDED because in chrono rotations are body-relative
797 B_overline_A.block(0, 3 + 6 * n, 3, 3) = block;
798
799 // delta epsilon_tilde_2_A
800 block = T_A[i].transpose() * L_alpha_beta_A[i](n, 1);
801 B_overline_A.block(3, 6 * n, 3, 3) = block;
802 block = T_A[i].transpose() * myA_2_X * Phi_Delta_A_n_LI_i;
803 block = block * this->m_nodes[n]->GetA(); //***NEEDED because in chrono rotations are body-relative
804 B_overline_A.block(3, 3 + 6 * n, 3, 3) = block;
805 }
806
807 B_overline_3_ABCD.row(i) = B_overline_A.row(2);
808 B_overline_6_ABCD.row(i) = B_overline_A.row(5);
809 }
810
811 for (int i = 0; i < NUMIP; i++) {
812 ChMatrixNM<double, 1, 4> sh1;
813 ChMatrixNM<double, 1, 4> sh2;
814 sh1(0, 0) = (1. + xi_i[i][1]) * 0.5;
815 sh1(0, 1) = 0;
816 sh1(0, 2) = (1. - xi_i[i][1]) * 0.5;
817 sh1(0, 3) = 0;
818
819 sh2(0, 0) = 0;
820 sh2(0, 1) = (1. - xi_i[i][0]) * 0.5;
821 sh2(0, 2) = 0;
822 sh2(0, 3) = (1. + xi_i[i][0]) * 0.5;
823
824 eps_tilde_1_i[i].z() = sh1(0, 0) * eps_tilde_1_A[0].z() + sh1(0, 2) * eps_tilde_1_A[2].z();
825 eps_tilde_2_i[i].z() = sh2(0, 1) * eps_tilde_2_A[1].z() + sh2(0, 3) * eps_tilde_2_A[3].z();
826
827 B_overline_i[i].row(2) = sh1 * B_overline_3_ABCD;
828 B_overline_i[i].row(5) = sh2 * B_overline_6_ABCD;
829 }
830
831 #endif
832
833 // EAS: B membranali
834 // {
835 // int tmpidx1[5] = {0, 1, 5, 4, 2};
836 // for (int i = 0; i < NUMIP; i++) {
837 // for (int n = 1; n <= 4; n++) {
838 //#if 0
839 // CopyMatrixRow(B_overline_m_i[i], n, B_overline_i[i], tmpidx1[n]);
840 //#endif
841 // B_overline_m_i[i].CopyMatrixRow(n, B_overline_i[i], tmpidx1[n]);
842 // }
843 // }
844 // }
845
846 /* Calcola le azioni interne */
847 for (int i = 0; i < NUMIP; i++) {
848 epsilon.segment(0, 3) = eps_tilde_1_i[i].eigen();
849 epsilon.segment(3, 3) = eps_tilde_2_i[i].eigen();
850 epsilon.segment(6, 3) = k_tilde_1_i[i].eigen();
851 epsilon.segment(9, 3) = k_tilde_2_i[i].eigen();
852
853 //***TODO*** add the EAS effect using the epsilon_hat
854 // epsilon_hat = P_i[i] * beta;
855 // epsilon += epsilon_hat;
856
857 ChVector<> eps_tot_1(epsilon.segment(0, 3));
858 ChVector<> eps_tot_2(epsilon.segment(3, 3));
859 ChVector<> k_tot_1(epsilon.segment(6, 3));
860 ChVector<> k_tot_2(epsilon.segment(9, 3));
861
862 // Compute strains using
863 // constitutive law of material
864
865 ChVector<> n1, n2, m1, m2;
866 ChVector<> l_n1, l_n2, l_m1, l_m2;
867 // loop on layers
868 for (size_t il = 0; il < this->m_layers.size(); ++il) {
869 // compute layer stresses (per-unit-length forces and torques), and accumulate
870 m_layers[il].GetMaterial()->ComputeStress(l_n1, l_n2, l_m1, l_m2, eps_tot_1, eps_tot_2, k_tot_1, k_tot_2,
871 m_layers_z[il], m_layers_z[il + 1], m_layers[il].Get_theta());
872 n1 += l_n1;
873 n2 += l_n2;
874 m1 += l_m1;
875 m2 += l_m2;
876 }
877
878 stress_i[i].segment(0,3) = n1.eigen();
879 stress_i[i].segment(3,3) = n2.eigen();
880 stress_i[i].segment(6,3) = m1.eigen();
881 stress_i[i].segment(9,3) = m2.eigen();
882
883 ChMatrix33<> Hh;
884 ChVector<> Tn1 = T_i[i] * n1;
885 ChVector<> Tn2 = T_i[i] * n2;
886 ChVector<> Tm1 = T_i[i] * m1;
887 ChVector<> Tm2 = T_i[i] * m2;
888 Hh = TensorProduct(Tn1, y_i_1[i]) - ChMatrix33<>(Tn1 ^ y_i_1[i]) + TensorProduct(Tn2, y_i_2[i]) -
889 ChMatrix33<>(Tn2 ^ y_i_2[i]) + TensorProduct(Tm1, k_1_i[i]) - ChMatrix33<>(Tm1 ^ k_1_i[i]) +
890 TensorProduct(Tm2, k_2_i[i]) - ChMatrix33<>(Tm2 ^ k_2_i[i]);
891
892 ChStarMatrix33<> Tn1_tilde(Tn1);
893 ChStarMatrix33<> Tn2_tilde(Tn2);
894 ChStarMatrix33<> Tm1_tilde(Tm1);
895 ChStarMatrix33<> Tm2_tilde(Tm2);
896
897 G_i[i].block(12, 0, 3, 3) = Tn1_tilde;
898 G_i[i].block(12, 3, 3, 3) = Tn2_tilde;
899 G_i[i].block(12, 6, 3, 3) = Tm1_tilde;
900 G_i[i].block(12, 9, 3, 3) = Tm2_tilde;
901 G_i[i].block(0, 12, 3, 3) = Tn1_tilde.transpose();
902 G_i[i].block(3, 12, 3, 3) = Tn2_tilde.transpose();
903 G_i[i].block(12, 12, 3, 3) = Hh;
904 }
905
906 // Residual
907
908 ChVectorN<double, 24> rd;
909 for (int i = 0; i < NUMIP; i++) {
910 rd = (-alpha_i[i] * w_i[i]) * B_overline_i[i].transpose() * stress_i[i];
911 Fi.segment(0, 24) += rd;
912
913 #ifdef CHUSE_EAS
914 double dCoef = 1.0; //***TODO*** autoset this
915 ChVectorN<double, IDOFS> rbeta;
916 rbeta = (-alpha_i[i] * w_i[i] / dCoef) * P_i[i].transpose() * stress_i[i];
917 Fi.segment(24, IDOFS) = rbeta;
918 #endif
919 }
920
921 // Add damping:
922
923 // fills velocities with speeds of nodes:
924 // {dxdt_1, angular_vel_1, dxdt_2, angular_vel_2, ...}
925 // Note that angular velocities are in node local frame because in C::E angular
926 // increments are assumed in local frame csys.
927 ChStateDelta velocities(24, nullptr);
928 this->LoadableGetStateBlock_w(0, velocities);
929
930 // Note that, in case of ChDampingReissnerRayleigh (rayleigh damping), this
931 // is like computing Fi_damping = -beta*[Km]*v without explicitly building [Km],
932 // rather exploit the factorization B'CB, as:
933 // Alpha*[Km]*v = beta * [sum_integraton_points ( [B_i]'*[C]*[B_i] * a_i * w_i )] *v
934 // Alpha*[Km]*v = sum_integraton_points ( [B_i]' * (beta * [C]*[B_i] *v) * a_i * w_i )
935 // Alpha*[Km]*v = sum_integraton_points ( [B_i]' * ( stress_damp_i ) * a_i * w_i )
936 // where for convenience we denoted stress_damp = beta*[C]*[B_i]*v
937
938 ChVectorN<double, 12> strain_dt;
939 ChVectorN<double, 12> stress_damp;
940
941 for (int i = 0; i < NUMIP; i++) {
942
943 strain_dt = B_overline_i[i] * velocities; // [B_i] *v
944
945 ChVector<> eps_dt_1(strain_dt.segment(0, 3));
946 ChVector<> eps_dt_2(strain_dt.segment(3, 3));
947 ChVector<> k_dt_1(strain_dt.segment(6, 3));
948 ChVector<> k_dt_2(strain_dt.segment(9, 3));
949 ChVector<> n1, n2, m1, m2;
950 ChVector<> l_n1, l_n2, l_m1, l_m2;
951 // loop on layers
952 for (size_t il = 0; il < this->m_layers.size(); ++il) {
953
954 if (m_layers[il].GetMaterial()->GetDamping()) {
955
956 // compute layer stresses (per-unit-length forces and torques), and accumulate [C]*[B_i]*v
957 m_layers[il].GetMaterial()->GetDamping()->ComputeStress(l_n1, l_n2, l_m1, l_m2, eps_dt_1, eps_dt_2, k_dt_1, k_dt_2,
958 m_layers_z[il], m_layers_z[il + 1], m_layers[il].Get_theta());
959 n1 += l_n1;
960 n2 += l_n2;
961 m1 += l_m1;
962 m2 += l_m2;
963 }
964 }
965 stress_damp.segment(0, 3) = n1.eigen();
966 stress_damp.segment(3, 3) = n2.eigen();
967 stress_damp.segment(6, 3) = m1.eigen();
968 stress_damp.segment(9, 3) = m2.eigen();
969
970 Fi.segment(0, 24) += (-alpha_i[i] * w_i[i]) * B_overline_i[i].transpose() * stress_damp;
971 }
972
973
974 }
975
976 // -----------------------------------------------------------------------------
977 // Jacobians of internal forces
978 // -----------------------------------------------------------------------------
979
ComputeInternalJacobians(double Kfactor,double Rfactor)980 void ChElementShellReissner4::ComputeInternalJacobians(double Kfactor, double Rfactor) {
981 m_JacobianMatrix.setZero();
982
983 // tangente
984
985 ChMatrixNM<double, 24, 24> Kg;
986 ChMatrixNM<double, 24, 24> Km;
987 ChMatrixNM<double, 24, 24> Rm;
988 ChMatrixNM<double, IDOFS, 24> K_beta_q;
989 ChMatrixNM<double, IDOFS, IDOFS> K_beta_beta;
990
991 for (int i = 0; i < NUMIP; i++) {
992 // GEOMETRIC STIFFNESS Kg:
993
994 Kg = D_overline_i[i].transpose() * G_i[i] * D_overline_i[i];
995
996 // MATERIAL STIFFNESS Km:
997
998 ChMatrixNM<double, 12, 12> C;
999 ChMatrixNM<double, 12, 12> l_C;
1000 C.setZero();
1001 l_C.setZero();
1002 // loop on layers
1003 for (size_t il = 0; il < m_layers.size(); ++il) {
1004 // compute layer tang. material stiff, and accumulate
1005 m_layers[il].GetMaterial()->ComputeStiffnessMatrix(
1006 l_C, eps_tilde_1_i[i], eps_tilde_2_i[i], k_tilde_1_i[i], k_tilde_2_i[i], m_layers_z[il],
1007 m_layers_z[il + 1],
1008 m_layers[il].Get_theta()); // ***TODO*** use the total epsilon including the 'hat' component from EAS
1009 C += l_C;
1010 }
1011
1012 Km = B_overline_i[i].transpose() * C * B_overline_i[i];
1013
1014 // EAS STIFFNESS K_beta_q:
1015 K_beta_q = P_i[i].transpose() * C * B_overline_i[i];
1016
1017 // EAS STIFFNESS K_beta_beta:
1018 K_beta_beta = P_i[i].transpose() * C * P_i[i];
1019
1020 // Assembly the entire jacobian
1021 // [ Km Kbq' ]
1022 // [ Kbq Kbb ]
1023
1024 double dCoef = 1.0; //***TODO*** autoset this
1025
1026 Km *= (alpha_i[i] * w_i[i] * dCoef * Kfactor); // was: ... * dCoef * (Kfactor+Rfactor * m_Alpha));
1027 m_JacobianMatrix += Km;
1028
1029 #ifdef CHUSE_KGEOMETRIC
1030 Kg *= (alpha_i[i] * w_i[i] * dCoef * Kfactor);
1031 m_JacobianMatrix += Kg;
1032 #endif
1033
1034 #ifdef CHUSE_EAS
1035 K_beta_q *= (alpha_i[i] * w_i[i] * Kfactor);
1036 K_beta_beta *= (alpha_i[i] * (w_i[i] / dCoef) * Kfactor);
1037
1038 m_JacobianMatrix.block(24, 0, 7, 24) += K_beta_q;
1039 m_JacobianMatrix.block(0, 24, 24, 7) += K_beta_q.transpose();
1040 m_JacobianMatrix.block(24, 4, 7, 7) += K_beta_beta;
1041 #endif
1042
1043 // Damping matrix
1044
1045 C.setZero();
1046 l_C.setZero();
1047 // loop on layers
1048 for (size_t il = 0; il < m_layers.size(); ++il) {
1049 // compute layer damping matrix
1050 if (m_layers[il].GetMaterial()->GetDamping()) {
1051 m_layers[il].GetMaterial()->GetDamping()->ComputeDampingMatrix(
1052 l_C,
1053 VNULL, VNULL, VNULL, VNULL, //***TODO*** should be more general: eps_dt_tilde_1_i[i], eps_dt_tilde_2_i[i], k_dt_tilde_1_i[i], k_dt_tilde_2_i[i],
1054 m_layers_z[il],
1055 m_layers_z[il + 1],
1056 m_layers[il].Get_theta());
1057 C += l_C;
1058 }
1059 }
1060 Rm = B_overline_i[i].transpose() * C * B_overline_i[i];
1061 Rm *= (alpha_i[i] * w_i[i] * dCoef * Rfactor);
1062 m_JacobianMatrix += Rm;
1063 }
1064 }
1065
1066 // -----------------------------------------------------------------------------
1067 // Shape functions
1068 // -----------------------------------------------------------------------------
1069
ShapeFunctions(ShapeVector & N,double x,double y)1070 void ChElementShellReissner4::ShapeFunctions(ShapeVector& N, double x, double y) {
1071 double xi[2];
1072 xi[0] = x, xi[1] = y;
1073
1074 N(0) = L1(xi);
1075 N(1) = L2(xi);
1076 N(2) = L3(xi);
1077 N(3) = L4(xi);
1078 }
1079
ShapeFunctionsDerivativeX(ShapeVector & Nx,double x,double y)1080 void ChElementShellReissner4::ShapeFunctionsDerivativeX(ShapeVector& Nx, double x, double y) {
1081 double xi[2];
1082 xi[0] = x, xi[1] = y;
1083 Nx(0) = L1_1(xi);
1084 Nx(1) = L2_1(xi);
1085 Nx(2) = L3_1(xi);
1086 Nx(3) = L4_1(xi);
1087 }
1088
ShapeFunctionsDerivativeY(ShapeVector & Ny,double x,double y)1089 void ChElementShellReissner4::ShapeFunctionsDerivativeY(ShapeVector& Ny, double x, double y) {
1090 double xi[2];
1091 xi[0] = x, xi[1] = y;
1092 Ny(0) = L1_2(xi);
1093 Ny(1) = L2_2(xi);
1094 Ny(2) = L3_2(xi);
1095 Ny(3) = L4_2(xi);
1096 }
1097
1098 // -----------------------------------------------------------------------------
1099 // Helper functions
1100 // -----------------------------------------------------------------------------
1101
1102 // -----------------------------------------------------------------------------
1103 // Interface to ChElementShell base class
1104 // -----------------------------------------------------------------------------
1105
EvaluateSectionDisplacement(const double u,const double v,ChVector<> & u_displ,ChVector<> & u_rotaz)1106 void ChElementShellReissner4::EvaluateSectionDisplacement(const double u,
1107 const double v,
1108 ChVector<>& u_displ,
1109 ChVector<>& u_rotaz) {
1110 // this is not a corotational element, so just do:
1111 EvaluateSectionPoint(u, v, u_displ);
1112 u_rotaz = VNULL; // no angles.. this is ANCF (or maybe return here the slope derivatives?)
1113 }
1114
EvaluateSectionFrame(const double u,const double v,ChVector<> & point,ChQuaternion<> & rot)1115 void ChElementShellReissner4::EvaluateSectionFrame(const double u,
1116 const double v,
1117 ChVector<>& point,
1118 ChQuaternion<>& rot) {
1119 // this is not a corotational element, so just do:
1120 EvaluateSectionPoint(u, v, point);
1121 rot = QUNIT; // or maybe use gram-schmidt to get csys of section from slopes?
1122 }
1123
EvaluateSectionPoint(const double u,const double v,ChVector<> & point)1124 void ChElementShellReissner4::EvaluateSectionPoint(const double u,
1125 const double v,
1126 ChVector<>& point) {
1127 ShapeVector N;
1128 this->ShapeFunctions(N, u, v);
1129
1130 const ChVector<>& pA = m_nodes[0]->GetPos();
1131 const ChVector<>& pB = m_nodes[1]->GetPos();
1132 const ChVector<>& pC = m_nodes[2]->GetPos();
1133 const ChVector<>& pD = m_nodes[3]->GetPos();
1134
1135 point = N(0) * pA + N(1) * pB + N(2) * pC + N(3) * pD;
1136 }
1137
1138 // -----------------------------------------------------------------------------
1139 // Functions for ChLoadable interface
1140 // -----------------------------------------------------------------------------
1141
1142 // Gets all the DOFs packed in a single vector (position part).
LoadableGetStateBlock_x(int block_offset,ChState & mD)1143 void ChElementShellReissner4::LoadableGetStateBlock_x(int block_offset, ChState& mD) {
1144 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos().eigen();
1145 mD.segment(block_offset + 3, 4) = m_nodes[0]->GetRot().eigen();
1146 mD.segment(block_offset + 7, 3) = m_nodes[1]->GetPos().eigen();
1147 mD.segment(block_offset + 10, 4) = m_nodes[1]->GetRot().eigen();
1148 mD.segment(block_offset + 14, 3) = m_nodes[2]->GetPos().eigen();
1149 mD.segment(block_offset + 17, 4) = m_nodes[2]->GetRot().eigen();
1150 mD.segment(block_offset + 21, 3) = m_nodes[3]->GetPos().eigen();
1151 mD.segment(block_offset + 24, 4) = m_nodes[3]->GetRot().eigen();
1152 }
1153
1154 // Gets all the DOFs packed in a single vector (velocity part).
LoadableGetStateBlock_w(int block_offset,ChStateDelta & mD)1155 void ChElementShellReissner4::LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) {
1156 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos_dt().eigen();
1157 mD.segment(block_offset + 3, 3) = m_nodes[0]->GetWvel_loc().eigen();
1158 mD.segment(block_offset + 6, 3) = m_nodes[1]->GetPos_dt().eigen();
1159 mD.segment(block_offset + 9, 3) = m_nodes[1]->GetWvel_loc().eigen();
1160 mD.segment(block_offset + 12, 3) = m_nodes[2]->GetPos_dt().eigen();
1161 mD.segment(block_offset + 15, 3) = m_nodes[2]->GetWvel_loc().eigen();
1162 mD.segment(block_offset + 18, 3) = m_nodes[3]->GetPos_dt().eigen();
1163 mD.segment(block_offset + 21, 3) = m_nodes[3]->GetWvel_loc().eigen();
1164 }
1165
LoadableStateIncrement(const unsigned int off_x,ChState & x_new,const ChState & x,const unsigned int off_v,const ChStateDelta & Dv)1166 void ChElementShellReissner4::LoadableStateIncrement(const unsigned int off_x, ChState& x_new, const ChState& x, const unsigned int off_v, const ChStateDelta& Dv) {
1167 for (int i = 0; i < 4; i++) {
1168 this->m_nodes[i]->NodeIntStateIncrement(off_x + 7 * i , x_new, x, off_v + 6 * i , Dv);
1169 }
1170 }
1171
EvaluateSectionVelNorm(double U,double V,ChVector<> & Result)1172 void ChElementShellReissner4::EvaluateSectionVelNorm(double U, double V, ChVector<>& Result) {
1173 ShapeVector N;
1174 ShapeFunctions(N, U, V);
1175 for (unsigned int ii = 0; ii < 4; ii++) {
1176 Result += N(ii) * m_nodes[ii]->GetPos_dt();
1177 }
1178 }
1179
1180 // Get the pointers to the contained ChVariables, appending to the mvars vector.
LoadableGetVariables(std::vector<ChVariables * > & mvars)1181 void ChElementShellReissner4::LoadableGetVariables(std::vector<ChVariables*>& mvars) {
1182 for (int i = 0; i < m_nodes.size(); ++i) {
1183 mvars.push_back(&m_nodes[i]->Variables());
1184 }
1185 }
1186
1187 // Evaluate N'*F , where N is the shape function evaluated at (U,V) coordinates of the surface.
ComputeNF(const double U,const double V,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1188 void ChElementShellReissner4::ComputeNF(
1189 const double U, // parametric coordinate in surface
1190 const double V, // parametric coordinate in surface
1191 ChVectorDynamic<>& Qi, // Return result of Q = N'*F here
1192 double& detJ, // Return det[J] here
1193 const ChVectorDynamic<>& F, // Input F vector, size is =n. field coords.
1194 ChVectorDynamic<>* state_x, // if != 0, update state (pos. part) to this, then evaluate Q
1195 ChVectorDynamic<>* state_w // if != 0, update state (speed part) to this, then evaluate Q
1196 ) {
1197 ShapeVector N;
1198 ShapeFunctions(N, U, V);
1199
1200 //***TODO*** exact det of jacobian at u,v
1201 detJ = GetLengthX() * GetLengthY() / 4.0; // approximate
1202
1203 Qi.segment(0, 3) = N(0) * F.segment(0, 3);
1204 Qi.segment(3, 3) = N(0) * F.segment(3, 3);
1205
1206 Qi.segment(6, 3) = N(1) * F.segment(0, 3);
1207 Qi.segment(9, 3) = N(1) * F.segment(3, 3);
1208
1209 Qi.segment(12, 3) = N(2) * F.segment(0, 3);
1210 Qi.segment(15, 3) = N(2) * F.segment(3, 3);
1211
1212 Qi.segment(18, 3) = N(3) * F.segment(0, 3);
1213 Qi.segment(21, 3) = N(3) * F.segment(3, 3);
1214 }
1215
1216 // Evaluate N'*F , where N is the shape function evaluated at (U,V,W) coordinates of the surface.
ComputeNF(const double U,const double V,const double W,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1217 void ChElementShellReissner4::ComputeNF(
1218 const double U, // parametric coordinate in volume
1219 const double V, // parametric coordinate in volume
1220 const double W, // parametric coordinate in volume
1221 ChVectorDynamic<>& Qi, // Return result of N'*F here, maybe with offset block_offset
1222 double& detJ, // Return det[J] here
1223 const ChVectorDynamic<>& F, // Input F vector, size is = n.field coords.
1224 ChVectorDynamic<>* state_x, // if != 0, update state (pos. part) to this, then evaluate Q
1225 ChVectorDynamic<>* state_w // if != 0, update state (speed part) to this, then evaluate Q
1226 ) {
1227 ShapeVector N;
1228 ShapeFunctions(N, U, V);
1229
1230 detJ = GetLengthX() * GetLengthY() / 4.0; // approximate
1231 detJ *= GetThickness();
1232
1233 Qi.segment(0,3) = N(0) * F.segment(0,3);
1234 Qi.segment(3,3) = N(0) * F.segment(3,3);
1235
1236 Qi.segment(6, 3) = N(1) * F.segment(0, 3);
1237 Qi.segment(9, 3) = N(1) * F.segment(3, 3);
1238
1239 Qi.segment(12, 3) = N(2) * F.segment(0, 3);
1240 Qi.segment(15, 3) = N(2) * F.segment(3, 3);
1241
1242 Qi.segment(18, 3) = N(3) * F.segment(0, 3);
1243 Qi.segment(21, 3) = N(3) * F.segment(3, 3);
1244 }
1245
1246 // -----------------------------------------------------------------------------
1247 //
1248 // -----------------------------------------------------------------------------
1249
1250 // Calculate average element density (needed for ChLoaderVolumeGravity).
GetDensity()1251 double ChElementShellReissner4::GetDensity() {
1252 double tot_density = 0;
1253 for (size_t kl = 0; kl < m_layers.size(); kl++) {
1254 double rho = m_layers[kl].GetMaterial()->Get_rho();
1255 double layerthick = m_layers[kl].Get_thickness();
1256 tot_density += rho * layerthick;
1257 }
1258 return tot_density / tot_thickness;
1259 }
1260
1261 // Calculate normal to the surface at (U,V) coordinates.
ComputeNormal(const double U,const double V)1262 ChVector<> ChElementShellReissner4::ComputeNormal(const double U, const double V) {
1263 ShapeVector Nx;
1264 ShapeVector Ny;
1265 ShapeFunctionsDerivativeX(Nx, U, V);
1266 ShapeFunctionsDerivativeY(Ny, U, V);
1267
1268 const ChVector<>& pA = m_nodes[0]->GetPos();
1269 const ChVector<>& pB = m_nodes[1]->GetPos();
1270 const ChVector<>& pC = m_nodes[2]->GetPos();
1271 const ChVector<>& pD = m_nodes[3]->GetPos();
1272
1273 ChVector<> Gx = Nx(0) * pA + Nx(1) * pB + Nx(2) * pC + Nx(3) * pD;
1274 ChVector<> Gy = Ny(0) * pA + Ny(1) * pB + Ny(2) * pC + Ny(3) * pD;
1275
1276 ChVector<> mnorm = Vcross(Gx, Gy);
1277 return mnorm.GetNormalized();
1278 }
1279
1280 // Private constructor (a layer can be created only by adding it to an element)
Layer(ChElementShellReissner4 * element,double thickness,double theta,std::shared_ptr<ChMaterialShellReissner> material)1281 ChElementShellReissner4::Layer::Layer(ChElementShellReissner4* element,
1282 double thickness,
1283 double theta,
1284 std::shared_ptr<ChMaterialShellReissner> material)
1285 : m_element(element), m_thickness(thickness), m_theta(theta), m_material(material) {}
1286
1287 // Initial setup for this layer:
SetupInitial()1288 void ChElementShellReissner4::Layer::SetupInitial() {}
1289
1290 } // end namespace fea
1291 } // end namespace chrono
1292