MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
shelleas.cc
Go to the documentation of this file.
1 /*
2  * MBDyn (C) is a multibody analysis code.
3  * http://www.mbdyn.org
4  *
5  * Copyright (C) 2010-2017
6  *
7  * Marco Morandini <morandini@aero.polimi.it>
8  * Riccardo Vescovini <vescovini@aero.polimi.it>
9  *
10  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11  * via La Masa, 34 - 20156 Milano, Italy
12  * http://www.aero.polimi.it
13  *
14  * Changing this copyright notice is forbidden.
15  *
16  * This program is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation (version 2 of the License).
19  *
20  *
21  * This program is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with this program; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  */
30 
31 /*
32  * Inspired by
33  * Wojciech Witkowski
34  * "4-Node combined shell element with semi-EAS-ANS strain interpolations
35  * in 6-parameter shell theories with drilling degrees of freedom"
36  * Comput Mech (2009) 43:307­319 DOI 10.1007/s00466-008-0307-x
37  */
38 
39 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
40 
41 #include "constltp_impl.h"
42 #include "tpldrive_impl.h"
43 #include "shelleas.h"
44 #include "mynewmem.h"
45 
46 #include "shell.hc"
47 
48 
49 /***************************************
50 *
51 *
52 ****************************************/
53 
54 doublereal Shell4EAS::xi_i[Shell4EAS::NUMIP][2] = {
55  {-1. / std::sqrt(3.), -1. / std::sqrt(3.)},
56  { 1. / std::sqrt(3.), -1. / std::sqrt(3.)},
57  { 1. / std::sqrt(3.), 1. / std::sqrt(3.)},
58  {-1. / std::sqrt(3.), 1. / std::sqrt(3.)}
59 
60 
61 // {-0.774596669241483, -0.774596669241483},
62 // {-0.774596669241483, 0.},
63 // {-0.774596669241483, 0.774596669241483},
64 // { 0., -0.774596669241483},
65 // { 0., 0.},
66 // { 0., 0.774596669241483},
67 // { 0.774596669241483, -0.774596669241483},
68 // { 0.774596669241483, 0.},
69 // { 0.774596669241483, 0.774596669241483}
70 
71 };
72 
73 doublereal Shell4EAS::w_i[Shell4EAS::NUMIP] = {
74  1.,
75  1.,
76  1.,
77  1.
78 
79 // 0.555555555555556 * 0.555555555555556,
80 // 0.555555555555556 * 0.888888888888889,
81 // 0.555555555555556 * 0.555555555555556,
82 // 0.888888888888889 * 0.555555555555556,
83 // 0.888888888888889 * 0.888888888888889,
84 // 0.888888888888889 * 0.555555555555556,
85 // 0.555555555555556 * 0.555555555555556,
86 // 0.555555555555556 * 0.888888888888889,
87 // 0.555555555555556 * 0.555555555555556
88 
89 };
90 
91 doublereal Shell4EAS::xi_n[Shell4EAS::NUMNODES][2] = {
92  { 1., 1.},
93  {-1., 1.},
94  {-1., -1.},
95  { 1., -1.}
96 };
97 
98 doublereal Shell4EAS::xi_0[2] = {0., 0.};
99 
100 void
101 Shell4EAS::UpdateNodalAndAveragePosAndOrientation(void)
102 {
103  Mat3x3 T_avg(Zero3x3);
104  Mat3x3 Tn[NUMNODES];
105  Mat3x3 R_tilde_n[NUMNODES];
106  for (integer i = 0; i < NUMNODES; i++) {
107 // xa[i] = pNode[i]->GetXCurr()+pNode[i]->GetRCurr()*f[i];
108  xa[i] = pNode[i]->GetXCurr();
109  Tn[i] = pNode[i]->GetRCurr() * iTa[i];
110 // T_a_avg += Ta[i];
111  T_avg += pNode[i]->GetRRef() * iTa[i];
112  }
113  T_avg /= 4.;
114  //FIXME; alternative solution: polar decomposition of T_a_avg = T0*U
115  T_overline = RotManip::Rot(RotManip::VecRot(T_avg));
116  for (integer i = 0; i < NUMNODES; i++) {
117  R_tilde_n[i] = T_overline.MulTM(Tn[i]);
118  phi_tilde_n[i] = RotManip::VecRot(R_tilde_n[i]);
119  }
120 }
121 
122 void
123 Shell4EAS::ComputeInitialNodeAndIptOrientation(void)
124 {
125  for (integer i = 0; i < NUMNODES; i++) {
126  xa[i] = pNode[i]->GetXCurr();
127  }
128  for (integer i = 0; i < NUMNODES; i++) {
129  Vec3 t1 = InterpDeriv_xi1(xa, xi_n[i]);
130  t1 = t1 / t1.Norm();
131  Vec3 t2 = InterpDeriv_xi2(xa, xi_n[i]);
132  t2 = t2 / t2.Norm();
133  Vec3 t3 = t1.Cross(t2);
134  t3 = t3 / t3.Norm();
135  t2 = t3.Cross(t1);
136  iTa[i] = (pNode[i]->GetRCurr()).MulTM(Mat3x3(t1, t2, t3));
137  }
138  for (integer i = 0; i < NUMIP; i++) {
139  iTa_i[i] = Eye3;
140  }
141  UpdateNodalAndAveragePosAndOrientation();
142  InterpolateOrientation();
143  for (integer i = 0; i < NUMIP; i++) {
144  Vec3 t1 = InterpDeriv_xi1(xa, xi_i[i]);
145  t1 = t1 / t1.Norm();
146  Vec3 t2 = InterpDeriv_xi2(xa, xi_i[i]);
147  t2 = t2 / t2.Norm();
148  Vec3 t3 = t1.Cross(t2);
149  t3 = t3 / t3.Norm();
150  t2 = t3.Cross(t1);
151  iTa_i[i] = (T_i[i]).MulTM(Mat3x3(t1, t2, t3));
152  }
153  InterpolateOrientation();
154 }
155 
156 void
157 Shell4EAS::InterpolateOrientation(void)
158 {
159  Mat3x3 DRot_I_phi_tilde_n_MT_T_overline[NUMNODES];
160  for (integer n = 0; n < NUMNODES; n++) {
161  DRot_I_phi_tilde_n_MT_T_overline[n] =
162  RotManip::DRot_I(phi_tilde_n[n]).MulMT(T_overline);
163  }
164  Mat3x3 Ri, Gammai;
165  for (integer i = 0; i < NUMIP; i++) {
166  phi_tilde_i[i] = Interp(phi_tilde_n, xi_i[i]);
167  RotManip::RotAndDRot(phi_tilde_i[i], Ri, Gammai);
168  T_i[i] = T_overline * Ri * iTa_i[i];
169  Mat3x3 T_overline_Gamma_tilde_i(T_overline * Gammai);
170  for (int n = 0; n < NUMNODES; n++) {
171  Phi_Delta_i[i][n] = T_overline_Gamma_tilde_i *
172  DRot_I_phi_tilde_n_MT_T_overline[n];
173  }
174  }
175  Vec3 phi_tilde_0 = Interp(phi_tilde_n, xi_0);
176  T_0 = T_overline * RotManip::Rot(phi_tilde_0);
177 }
178 
179 // void
180 // Shell4EAS::ComputeIPSEPRotations(void)
181 // {
182 // for (integer i = 0; i < NUMIP; i++) {
183 // Q_i[i] = T_i[i].MulMT(T_0_i[i]);
184 // }
185 // for (integer i = 0; i < NUMSSEP; i++) {
186 // Q_A[i] = T_A[i].MulMT(T_0_A[i]);
187 // }
188 // }
189 
190 void
191 Shell4EAS::ComputeIPCurvature(void)
192 {
193  Mat3x3 Gamma_I_n_MT_T_overline[NUMNODES];
194  for (integer n = 0; n < NUMNODES; n++) {
195  Gamma_I_n_MT_T_overline[n] = RotManip::DRot_I(phi_tilde_n[n]).MulMT(T_overline);
196  }
197  for (integer i = 0; i < NUMIP; i++) {
198  Vec3 phi_tilde_1_i(Zero3);
199  Vec3 phi_tilde_2_i(Zero3);
200  InterpDeriv(phi_tilde_n, L_alpha_beta_i[i], phi_tilde_1_i, phi_tilde_2_i);
201  Mat3x3 T_overlineGamma_tilde_i(T_overline * RotManip::DRot(phi_tilde_i[i]));
202  k_1_i[i] = T_overlineGamma_tilde_i * phi_tilde_1_i;
203  k_2_i[i] = T_overlineGamma_tilde_i * phi_tilde_2_i;
204  Mat3x3 tmp1 = T_overline * RotManip::Elle(phi_tilde_i[i], phi_tilde_1_i);
205  Mat3x3 tmp2 = T_overline * RotManip::Elle(phi_tilde_i[i], phi_tilde_2_i);
206  for (int n = 0; n < NUMNODES; n++) {
207  Kappa_delta_i_1[i][n] = tmp1 * Gamma_I_n_MT_T_overline[n] * LI[n](xi_i[i]) +
208  Phi_Delta_i[i][n] * L_alpha_beta_i[i](n + 1, 1);
209  Kappa_delta_i_2[i][n] = tmp2 * Gamma_I_n_MT_T_overline[n] * LI[n](xi_i[i]) +
210  Phi_Delta_i[i][n] * L_alpha_beta_i[i](n + 1, 2);
211  }
212  }
213 }
214 
215 Shell4EAS::Shell4EAS(unsigned int uL,
216  const DofOwner* pDO,
217  const StructNode* pN1, const StructNode* pN2,
218  const StructNode* pN3, const StructNode* pN4,
219 #if 0 // TODO: offset
220  const Vec3& f1, const Vec3& f2,
221  const Vec3& f3, const Vec3& f4,
222 #endif
223  const Mat3x3& R1, const Mat3x3& R2,
224  const Mat3x3& R3, const Mat3x3& R4,
225 #ifdef USE_CL_IN_SHELL
226  const ConstitutiveLaw<vh, fmh>** pDTmp,
227 #else // ! USE_CL_IN_SHELL
228  const fmh& pDTmp,
229  const vh& PreStress,
230 #endif // ! USE_CL_IN_SHELL
231  flag fOut)
232 :
233 Elem(uL, fOut),
234 Shell(uL, pDO, fOut),
235 
236 S_alpha_beta_0(2, 2),
237 
238 L_alpha_beta_i(NUMIP, fmh(4, 2) ),
239 
240 B_overline_i(NUMIP, fmh(12, 24) ),
241 
242 P_i(NUMIP, fmh(12, iGetNumDof()) ),
243 
244 beta(iGetNumDof()),
245 epsilon_hat(12),
246 epsilon(12),
247 
248 #ifndef USE_CL_IN_SHELL
249 bPreStress(PreStress.Norm() > 0.),
250 PreStress(PreStress),
251 #endif // ! USE_CL_IN_SHELL
252 
253 DRef(NUMIP, fmh(12, 12)),
254 stress_i(NUMIP, vh(12))
255 {
256 #ifdef USE_CL_IN_SHELL
257  for (integer i = 0; i < NUMIP; i++) {
258  pD[i] = 0;
259  SAFENEWWITHCONSTRUCTOR(pD[i],
260  ConstitutiveLawOwnerType,
261  ConstitutiveLawOwnerType(pDTmp[i]));
262  }
263 #else // ! USE_CL_IN_SHELL
264  for (unsigned i = 0; i < NUMIP; i++) {
265  DRef[i] = pDTmp;
266  }
267 #endif // ! USE_CL_IN_SHELL
268 
269  pNode[NODE1] = pN1;
270  pNode[NODE2] = pN2;
271  pNode[NODE3] = pN3;
272  pNode[NODE4] = pN4;
273 // f[NODE1] = f1;
274 // f[NODE2] = f2;
275 // f[NODE3] = f3;
276 // f[NODE4] = f4;
277  ComputeInitialNodeAndIptOrientation();
278 // UpdateNodalAndAveragePosAndOrientation();
279 // InterpolateOrientation();
280  // copy ref values
281  T0_overline = T_overline;
282  T_0_0 = T_0;
283  for (integer i = 0; i < NUMNODES; i++) {
284  xa_0[i] = xa[i];
285  }
286  for (integer i = 0; i < NUMIP; i++) {
287  T_0_i[i] = T_i[i];
288  }
289 
290 
291  fmh M_0(6, 6);
292  fmh M_0_Inv(6, 6);
293  {
294  Vec3 x_1 = InterpDeriv_xi1(xa, xi_0);
295  Vec3 x_2 = InterpDeriv_xi2(xa, xi_0);
296  S_alpha_beta_0(1, 1) = T_0_0.GetCol(1) * x_1;
297  S_alpha_beta_0(2, 1) = T_0_0.GetCol(2) * x_1;
298  S_alpha_beta_0(1, 2) = T_0_0.GetCol(1) * x_2;
299  S_alpha_beta_0(2, 2) = T_0_0.GetCol(2) * x_2;
300  alpha_0 = S_alpha_beta_0(1, 1) * S_alpha_beta_0(2, 2) -
301  S_alpha_beta_0(1, 2) * S_alpha_beta_0(2, 1);
302 
303  M_0(1, 1) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(1, 1);
304  M_0(1, 2) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(1, 2);
305  M_0(1, 3) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(1, 1);
306  M_0(1, 4) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(1, 2);
307 
308  M_0(2, 1) = S_alpha_beta_0(2, 1) * S_alpha_beta_0(2, 1);
309  M_0(2, 2) = S_alpha_beta_0(2, 2) * S_alpha_beta_0(2, 2);
310  M_0(2, 3) = S_alpha_beta_0(2, 2) * S_alpha_beta_0(2, 1);
311  M_0(2, 4) = S_alpha_beta_0(2, 1) * S_alpha_beta_0(2, 2);
312 
313  M_0(3, 1) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(2, 1);
314  M_0(3, 2) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(2, 2);
315  M_0(3, 3) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(2, 2);
316  M_0(3, 4) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(2, 1);
317 
318  M_0(4, 1) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(2, 1);
319  M_0(4, 2) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(2, 2);
320  M_0(4, 3) = S_alpha_beta_0(1, 2) * S_alpha_beta_0(2, 1);
321  M_0(4, 4) = S_alpha_beta_0(1, 1) * S_alpha_beta_0(2, 2);
322 
323  M_0(5, 5) = S_alpha_beta_0(1, 1);
324  M_0(5, 6) = S_alpha_beta_0(1, 2);
325  M_0(6, 5) = S_alpha_beta_0(2, 1);
326  M_0(6, 6) = S_alpha_beta_0(2, 2);
327 
328  InvBlockDiagonal4_2x4_2(M_0, M_0_Inv);
329  }
330 
331  for (integer i = 0; i < NUMIP; i++) {
332  fmh L_alpha_B_i(4, 2);
333  fmh S_alpha_beta_i(2, 2);
334  Vec3 x_1 = InterpDeriv_xi1(xa, xi_i[i]);
335  Vec3 x_2 = InterpDeriv_xi2(xa, xi_i[i]);
336  S_alpha_beta_i(1, 1) = T_0_i[i].GetCol(1) * x_1;
337  S_alpha_beta_i(2, 1) = T_0_i[i].GetCol(2) * x_1;
338  S_alpha_beta_i(1, 2) = T_0_i[i].GetCol(1) * x_2;
339  S_alpha_beta_i(2, 2) = T_0_i[i].GetCol(2) * x_2;
340  // alpha_i = det(S_alpha_beta_i)
341  alpha_i[i] = S_alpha_beta_i(1, 1) * S_alpha_beta_i(2, 2) -
342  S_alpha_beta_i(1, 2) * S_alpha_beta_i(2, 1);
343 
344  // xi_i_i = S_alpha_beta_i^{-1}
345  FullMatrixHandler xi_i_i(2, 2);
346  Inv2x2(S_alpha_beta_i, xi_i_i);
347 
348  for (integer n = 0; n < NUMNODES; n++) {
349  for (integer ii = 0; ii < 2; ii++) {
350  L_alpha_B_i(n + 1, ii + 1) = LI_J[n][ii](xi_i[i]);
351  }
352  }
353 
354  L_alpha_B_i.MatMatMul(L_alpha_beta_i[i], xi_i_i);
355 
356  fmh H(6, iGetNumDof());
357  doublereal t = xi_i[i][0] * xi_i[i][1];
358  H(1, 1) = xi_i[i][0];
359  H(1, 2) = t;
360 
361  H(2, 3) = xi_i[i][1];
362  H(2, 4) = t;
363 
364  H(3, 5) = xi_i[i][0];
365  H(3, 6) = t;
366 
367  H(4, 7) = xi_i[i][1];
368  H(4, 6) = t;
369 
370  doublereal a = 2. * xi_i[i][0] * (xi_i[i][1] * xi_i[i][1] - 1.);
371  doublereal b = 2. * xi_i[i][1] * (xi_i[i][0] * xi_i[i][0] - 1.);
372  H(5, 8) = a;
373  H(5, 9) = b;
374  H(5, 10) = a * b;
375  H(6, 11) = b;
376  H(6, 12) = a;
377  H(6, 13) = a * b;
378 // doublereal a = 2. * xi_i[i][0] * (xi_i[i][1] * xi_i[i][1] - 1.);
379 // doublereal b = 2. * xi_i[i][1] * (xi_i[i][0] * xi_i[i][0] - 1.);
380 // H(5, 9) = a;
381 // H(5, 10) = b;
382 // H(5, 11) = a * b;
383 // H(6, 12) = b;
384 // H(6, 13) = a;
385 // H(6, 14) = a * b;
386  fmh Perm(12, 6), tmpP(6, iGetNumDof());
387  // 1, 5, 4, 2, 3, 6
388  Perm(1, 1) = 1.;
389  Perm(2, 3) = 1.;
390  Perm(3, 5) = 1.;
391  Perm(4, 4) = 1.;
392  Perm(5, 2) = 1.;
393  Perm(6, 6) = 1.;
394 
395  M_0_Inv.MatTMatMul(tmpP, H);
396  Perm.MatMatMul(P_i[i], tmpP);
397  P_i[i].ScalarMul(alpha_0 / alpha_i[i]);
398  }
399  // save initial axial values
400  ComputeIPCurvature();
401  for (integer i = 0; i < NUMIP; i++) {
402  InterpDeriv(xa, L_alpha_beta_i[i], y_i_1[i], y_i_2[i]);
403  eps_tilde_1_0_i[i] = T_i[i].MulTV(y_i_1[i]);
404  eps_tilde_2_0_i[i] = T_i[i].MulTV(y_i_2[i]);
405  k_tilde_1_0_i[i] = T_i[i].MulTV(k_1_i[i]);
406  k_tilde_2_0_i[i] = T_i[i].MulTV(k_2_i[i]);
407  }
408 }
409 
410 Shell4EAS::~Shell4EAS(void)
411 {
412 #ifdef USE_CL_IN_SHELL
413  for (integer i = 0; i < NUMIP; i++) {
414  ASSERT(pD[i] != 0);
415  if (pD[i] != 0) {
416  SAFEDELETE(pD[i]);
417  }
418  }
419 #endif // USE_CL_IN_SHELL
420 }
421 
422 SubVectorHandler& Shell4EAS::AssRes(SubVectorHandler& WorkVec,
423  doublereal dCoef,
424  const VectorHandler& XCurr,
425  const VectorHandler& XPrimeCurr)
426 {
427  /* Dimensiona e resetta la matrice di lavoro */
428  integer iNumRows = 0;
429  integer iNumCols = 0;
430  WorkSpaceDim(&iNumRows, &iNumCols);
431  WorkVec.ResizeReset(iNumRows);
432 
433  /* Recupera e setta gli indici */
434  for (integer i = 0; i < 4; i++) {
435  integer iNodeFirstMomIndex = pNode[i]->iGetFirstMomentumIndex();
436  for (int iCnt = 1; iCnt <= 6; iCnt++) {
437  WorkVec.PutRowIndex(iCnt + 6 * i, iNodeFirstMomIndex + iCnt);
438  }
439  }
440 
441  integer iFirstReactionIndex = iGetFirstIndex();
442  for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
443  WorkVec.PutRowIndex(24 + iCnt, iFirstReactionIndex + iCnt);
444  }
445 
446  UpdateNodalAndAveragePosAndOrientation();
447  InterpolateOrientation();
448 // ComputeIPSEPRotations();
449  for (unsigned int i = 1; i <= iGetNumDof(); i++) {
450  beta(i) = XCurr(iFirstReactionIndex + i);
451  }
452 
453 
454  ComputeIPCurvature();
455  for (integer i = 0; i < NUMIP; i++) {
456  InterpDeriv(xa, L_alpha_beta_i[i], y_i_1[i], y_i_2[i]);
457  eps_tilde_1_i[i] = T_i[i].MulTV(y_i_1[i]) - eps_tilde_1_0_i[i];
458  eps_tilde_2_i[i] = T_i[i].MulTV(y_i_2[i]) - eps_tilde_2_0_i[i];
459  k_tilde_1_i[i] = T_i[i].MulTV(k_1_i[i]) - k_tilde_1_0_i[i];
460  k_tilde_2_i[i] = T_i[i].MulTV(k_2_i[i]) - k_tilde_2_0_i[i];
461  // parte variabile di B_overline_i
462  for (integer n = 0; n < NUMNODES; n++) {
463  Mat3x3 Phi_Delta_i_n_LI_i = Phi_Delta_i[i][n] * LI[n](xi_i[i]);
464 
465  // delta epsilon_tilde_1_i
466  B_overline_i[i].PutT(1, 1 + 6 * n, T_i[i] * L_alpha_beta_i[i](n + 1, 1));
467  B_overline_i[i].Put(1, 4 + 6 * n,
468  T_i[i].MulTM(Mat3x3(MatCross, y_i_1[i])) * Phi_Delta_i_n_LI_i
469  );
470 
471  // delta epsilon_tilde_2_i
472  B_overline_i[i].PutT(4, 1 + 6 * n, T_i[i] * L_alpha_beta_i[i](n + 1, 2));
473  B_overline_i[i].Put(4, 4 + 6 * n,
474  T_i[i].MulTM(Mat3x3(MatCross, y_i_2[i])) * Phi_Delta_i_n_LI_i
475  );
476 
477  // delta k_tilde_1_i
478  Vec3 phi_tilde_1_i(Zero3);
479  Vec3 phi_tilde_2_i(Zero3);
480  InterpDeriv(phi_tilde_n, L_alpha_beta_i[i], phi_tilde_1_i, phi_tilde_2_i);
481  B_overline_i[i].Put(7, 4 + 6 * n,
482  T_i[i].MulTM(Mat3x3(MatCross, k_1_i[i])) * Phi_Delta_i_n_LI_i
483  +
484  T_i[i].MulTM(Kappa_delta_i_1[i][n])
485  );
486 
487  // delta k_tilde_2_i
488  B_overline_i[i].Put(10, 4 + 6 * n,
489  T_i[i].MulTM(Mat3x3(MatCross, k_2_i[i])) * Phi_Delta_i_n_LI_i
490  +
491  T_i[i].MulTM(Kappa_delta_i_2[i][n])
492  );
493  }
494  }
495 
496 
497  /* Calcola le azioni interne */
498  for (integer i = 0; i < NUMIP; i++) {
499  epsilon.Put(1, eps_tilde_1_i[i]);
500  epsilon.Put(4, eps_tilde_2_i[i]);
501  epsilon.Put(7, k_tilde_1_i[i]);
502  epsilon.Put(10, k_tilde_2_i[i]);
503  // TODO: recupera epsilon_hat con l'ordine giusto per qua
504  P_i[i].MatVecMul(epsilon_hat, beta);
505 // {
506 // for (int off = 0; off < 4; off++) {
507 // Vec3 e1(epsilon_hat(1 + 3 * off), epsilon_hat(2 + 3 * off), epsilon_hat(3 + 3 * off));
508 // e1 = T_i[i].MulTV(e1);
509 // epsilon_hat(1 + 3 * off) = e1(1);
510 // epsilon_hat(2 + 3 * off) = e1(2);
511 // epsilon_hat(3 + 3 * off) = e1(3);
512 // }
513 // }
514 
515  epsilon += epsilon_hat;
516 #ifdef USE_CL_IN_SHELL
517  pD[i]->Update(epsilon);
518  stress_i[i] = pD[i]->GetF();
519 #else // ! USE_CL_IN_SHELL
520  DRef[i].MatVecMul(stress_i[i], epsilon);
521  if (bPreStress) {
522  stress_i[i] += PreStress;
523  }
524 #endif // ! USE_CL_IN_SHELL
525  }
526 
527 
528  //Residuo
529  //forze di volume
530  MyVectorHandler rd(24);
531  MyVectorHandler rbeta(iGetNumDof());
532  for (integer i = 0; i < NUMIP; i++) {
533  B_overline_i[i].MatTVecMul(rd, stress_i[i]);
534  P_i[i].MatTVecMul(rbeta, stress_i[i]);
535 
536  AssembleVector(WorkVec, 1, rd, -alpha_i[i] * w_i[i]);
537  AssembleVector(WorkVec, 25, rbeta, -alpha_i[i] * w_i[i]);
538  }
539 
540  return WorkVec;
541 }
542 #include <iostream>
543 // Jacobian matrix assembly
544 VariableSubMatrixHandler&
545 Shell4EAS::AssJac(VariableSubMatrixHandler& WorkMat,
546  doublereal dCoef,
547  const VectorHandler& XCurr,
548  const VectorHandler& XPrimeCurr)
549 {
550  FullSubMatrixHandler& WM = WorkMat.SetFull();
551 
552  /* Dimensiona e resetta la matrice di lavoro */
553  integer iNumRows = 0;
554  integer iNumCols = 0;
555  WorkSpaceDim(&iNumRows, &iNumCols);
556  WM.ResizeReset(iNumRows, iNumCols);
557 
558  /* Recupera e setta gli indici */
559  integer iFirstReactionIndex = iGetFirstIndex();
560  for (integer i = 0; i < 4; i++) {
561  integer iNodeFirstMomIndex = pNode[i]->iGetFirstMomentumIndex();
562  integer iNodeFirstPosIndex = pNode[i]->iGetFirstPositionIndex();
563  for (int iCnt = 1; iCnt <= 6; iCnt++) {
564  WM.PutRowIndex(iCnt + 6 * i, iNodeFirstMomIndex + iCnt);
565  WM.PutColIndex(iCnt + 6 * i, iNodeFirstPosIndex + iCnt);
566  }
567  }
568  for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
569  WM.PutRowIndex(24 + iCnt, iFirstReactionIndex + iCnt);
570  WM.PutColIndex(24 + iCnt, iFirstReactionIndex + iCnt);
571  }
572 
573 
578 // for (integer ip = 0; ip < NUMIP; ip++) {
579 // /* Finite difference check of B_overline matrix */
580 // Vec3 eps1, eps2, k1, k2, y1, y2;
581 //
582 // fmh B_overline_base(12, 24);
583 // B_overline_base = B_overline_i[ip];
584 // fmh pippo(12, 24);
585 // fmh pippodiff(12, 24);
586 //
587 // MyVectorHandler basedef(12);
588 // MyVectorHandler incdef(12);
589 //
590 // MyVectorHandler inc(50);
591 // MyVectorHandler xprime(50);
592 // inc.Reset();
593 // xprime.Reset();
594 //
595 // std::cerr << "----------- B_overline[" << ip << "]-----------" << std::endl;
596 // doublereal ddd = 0.000001;
597 // doublereal maxdiff = 0.;
598 // std::cerr << "\nxxxxxxxxxxxxxxx\n" << std::endl;
599 // std::cerr << std::fixed << std::setprecision(6) << B_overline_i[ip] << std::endl;
600 // for (integer n = 0; n < NUMNODES; n++) {
601 // Vec3 xbase = xa[n];
602 // Mat3x3 RDelta = pNode[n]->GetRCurr().MulMT(pNode[n]->GetRRef());
603 // Vec3 gbase(CGR_Rot::Param, RDelta);
604 // for (integer gdl = 1; gdl <= 3; gdl ++) {
605 // inc.PutCoef(pNode[n]->iGetFirstIndex() + 0 + gdl, xbase(gdl));
606 // inc.PutCoef(pNode[n]->iGetFirstIndex() + 3 + gdl, gbase(gdl));
607 // }
608 // }
609 //
610 // for (integer n = 0; n < NUMNODES; n++) {
611 // for (integer gdl = 1; gdl <= 6; gdl ++) {
612 // if (n == 1 && gdl == 1) {
613 // std::cerr << "Eccomi!\n";
614 // for (integer nn=0; nn < 4; nn++) {
615 // std::cerr << std::setprecision(12) << "node " << nn
616 // << " " << xa[nn] << "\n";
617 // }
618 // std::cerr << std::setprecision(12) << "etilde 1 " << eps_tilde_1_i[ip] << "\n";
619 // std::cerr << std::setprecision(12) << "etilde 2 " << eps_tilde_2_i[ip] << "\n";
620 // std::cerr << std::setprecision(12) << "y_i_1 " << y_i_1[ip] << "\n";
621 // std::cerr << std::setprecision(12) << "y_i_2 " << y_i_2[ip] << "\n";
622 // std::cerr << std::setprecision(12) << "T y_i_1 " << T_i[ip].MulTV(y_i_1[ip]) << "\n";
623 // std::cerr << std::setprecision(12) << "T y_i_2 " << T_i[ip].MulTV(y_i_2[ip]) << "\n";
624 // std::cerr << std::setprecision(12) << "T " << T_i[ip] << "\n";
625 // }
626 // doublereal saved_value = inc(pNode[n]->iGetFirstIndex() + gdl);
627 // inc.IncCoef(pNode[n]->iGetFirstIndex() + gdl, ddd);
628 // const_cast<StructNode*>(pNode[n])->Update(inc, xprime);
629 // UpdateNodalAndAveragePosAndOrientation();
630 // InterpolateOrientation();
631 // ComputeIPCurvature();
632 // InterpDeriv(xa, L_alpha_beta_i[ip], y1, y2);
633 // eps1 = T_i[ip].MulTV(y1) - Vec3(1., 0., 0.);
634 // eps2 = T_i[ip].MulTV(y2) - Vec3(0., 1., 0.);
635 // k1 = T_i[ip].MulTV(k_1_i[ip]) - T_0_i[ip].MulTV(k_1_0_i[ip]);
636 // k2 = T_i[ip].MulTV(k_2_i[ip]) - T_0_i[ip].MulTV(k_2_0_i[ip]);
637 // if (n == 1 && gdl == 1) {
638 // std::cerr << std::setprecision(12) << "e1 " << eps1 << "\n";
639 // std::cerr << std::setprecision(12) << "e2 " << eps2 << "\n";
640 // std::cerr << std::setprecision(12) << "y1 " << y1 << "\n";
641 // std::cerr << std::setprecision(12) << "y2 " << y2 << "\n";
642 // std::cerr << std::setprecision(12) << "T y_1 " << T_i[ip].MulTV(y1) << "\n";
643 // std::cerr << std::setprecision(12) << "T y_2 " << T_i[ip].MulTV(y2) << "\n";
644 // std::cerr << std::setprecision(12) << T_i[ip] << "\n";
645 // std::cerr << ",,,,,,,,,,,,,,,,,\n";
646 //
647 // }
648 // // Backup changes
649 // inc.PutCoef(pNode[n]->iGetFirstPositionIndex() + gdl, saved_value);
650 // const_cast<StructNode*>(pNode[n])->Update(inc, xprime);
651 // UpdateNodalAndAveragePosAndOrientation();
652 // InterpolateOrientation();
653 // ComputeIPCurvature();
654 // InterpDeriv(xa, L_alpha_beta_i[ip], y1, y2);
655 // for (integer comp = 1; comp <= 3; comp++) {
656 // doublereal d1 = (eps1(comp) - eps_tilde_1_i[ip](comp)) / ddd;
657 // doublereal d2 = (eps2(comp) - eps_tilde_2_i[ip](comp)) / ddd;
658 // doublereal d3 = (k1(comp) - k_tilde_1_i[ip](comp)) / ddd;
659 // doublereal d4 = (k2(comp) - k_tilde_2_i[ip](comp)) / ddd;
660 // d1 = std::abs(d1) > 1.E-100 ? d1 : 0.;
661 // d2 = std::abs(d2) > 1.E-100 ? d2 : 0.;
662 // d3 = std::abs(d3) > 1.E-100 ? d3 : 0.;
663 // d4 = std::abs(d4) > 1.E-100 ? d4 : 0.;
664 // pippo( comp, gdl + 6 * n) = d1;
665 // pippo(3 + comp, gdl + 6 * n) = d2;
666 // pippo(6 + comp, gdl + 6 * n) = d3;
667 // pippo(9 + comp, gdl + 6 * n) = d4;
668 //
669 // pippodiff(comp, gdl + 6 * n) =
670 // std::abs(pippo(comp, gdl + 6 * n) -
671 // B_overline_i[ip](comp, gdl + 6 * n)) > 1.E-30?
672 // pippo(comp, gdl + 6 * n) -
673 // B_overline_i[ip](comp, gdl + 6 * n) : 0.;
674 // pippodiff(3 + comp, gdl + 6 * n) =
675 // std::abs(pippo(3 + comp, gdl + 6 * n) -
676 // B_overline_i[ip](3 + comp, gdl + 6 * n) ) > 1.E-30?
677 // pippo(3 + comp, gdl + 6 * n) -
678 // B_overline_i[ip](3 + comp, gdl + 6 * n) : 0.;
679 // pippodiff(6 + comp, gdl + 6 * n) =
680 // std::abs(pippo(6 + comp, gdl + 6 * n) -
681 // B_overline_i[ip](6 + comp, gdl + 6 * n)) > 1.E-30?
682 // pippo(6 + comp, gdl + 6 * n) -
683 // B_overline_i[ip](6 + comp, gdl + 6 * n) : 0.;
684 // pippodiff(9 + comp, gdl + 6 * n) =
685 // std::abs(pippo(9 + comp, gdl + 6 * n) -
686 // B_overline_i[ip](9 + comp, gdl + 6 * n)) > 1.E-30?
687 // pippo(9 + comp, gdl + 6 * n) -
688 // B_overline_i[ip](9 + comp, gdl + 6 * n) : 0.;
689 // maxdiff = std::max(maxdiff, std::abs(pippodiff(0 + comp, gdl + 6 * n)));
690 // maxdiff = std::max(maxdiff, std::abs(pippodiff(3 + comp, gdl + 6 * n)));
691 // maxdiff = std::max(maxdiff, std::abs(pippodiff(6 + comp, gdl + 6 * n)));
692 // maxdiff = std::max(maxdiff, std::abs(pippodiff(9 + comp, gdl + 6 * n)));
693 // }
694 // }
695 // // check rotation linearization
696 // for (integer gdl = 4; gdl <= 6; gdl ++) {
697 // //Mat3x3 RDelta = pNode[n]->GetRCurr().MulMT(pNode[n]->GetRRef());
698 // //Mat3x3 Gamma = RotManip::DRot(RotManip::VecRot(RDelta));
699 // Mat3x3 Gamma = Mat3x3(CGR_Rot::MatG, pNode[n]->GetgCurr());
700 // doublereal saved_value = inc(pNode[n]->iGetFirstIndex() + gdl);
701 // Mat3x3 SavedT = T_i[ip];
702 // inc.IncCoef(pNode[n]->iGetFirstIndex() + gdl, ddd);
703 // const_cast<StructNode*>(pNode[n])->Update(inc, xprime);
704 // UpdateNodalAndAveragePosAndOrientation();
705 // InterpolateOrientation();
706 // ComputeIPCurvature();
707 // Vec3 phi_delta = ((T_i[ip]-SavedT).MulMT(SavedT)).Ax() / ddd;
708 // // Backup changes
709 // inc.PutCoef(pNode[n]->iGetFirstPositionIndex() + gdl, saved_value);
710 // const_cast<StructNode*>(pNode[n])->Update(inc, xprime);
711 // UpdateNodalAndAveragePosAndOrientation();
712 // InterpolateOrientation();
713 // ComputeIPCurvature();
714 // InterpDeriv(xa, L_alpha_beta_i[ip], y1, y2);
715 // std::cerr << std::fixed << std::setprecision(12) << "gdl "
716 // << gdl << " " << phi_delta << " ; "
717 // << (Phi_Delta_i[ip][n] * Gamma ).GetCol(gdl - 3) * LI[n](xi_i[ip]) << "\n";
718 // }
719 // }
720 // // for (integer j = 1; j <= pJac->iGetNumCols(); j++) {
721 // // pippo.PutCoef(j, i,
722 // // std::abs(incsol(j)) >
723 // // 1.E-100?incsol(j):0.);
724 // // pippodiff.PutCoef(j, i,
725 // // std::abs(pippo(j,i)-pJac->operator()(j,i)) >
726 // // 1.E-30?pippo(j,i)-pJac->operator()(j,i):0.);
727 // // maxdiff = std::max(maxdiff, std::abs(pippodiff(j, i)));
728 // // // use finite difference Jacobian matrix
729 // // pJac->PutCoef(j, i, incsol(j));
730 // // }
731 // std::cerr << "\n---------------\n" << std::endl;
732 // std::cerr << std::fixed << std::setprecision(6) << pippo << std::endl;
733 // std::cerr << "\n===============\n" << std::endl;
734 // std::cerr << std::fixed << std::setprecision(6) << pippodiff << std::endl;
735 // std::cerr << "\nZZZZZZZZZZZZZZZ\n" << std::endl;
736 // std::cerr << "MAXDIFF: " << std::fixed << std::setprecision(12) << maxdiff << std::endl;
737 // std::cerr << "\nXXXXXXXXXXXXXXX\n" << std::endl;
738 //
739 // }
744 
745 
746 
747 
748 
749 
750  // tangente
751 
752  FullMatrixHandler Km(24, 24);
753  FullMatrixHandler K_beta_q(iGetNumDof(), 24);
754  FullMatrixHandler K_beta_beta(iGetNumDof(), iGetNumDof());
755 
756  FullMatrixHandler Ktm(24, 12);
757 // FullMatrixHandler Kttm(24, 24);
758 
759  FullMatrixHandler Ktbetaq(iGetNumDof(), 12);
760 
761  FullMatrixHandler C(12, 12);
762  for (integer i = 0; i < NUMIP; i++) {
763 
764 #ifdef USE_CL_IN_SHELL
765  C = pD[i]->GetFDE();
766 #else // ! USE_CL_IN_SHELL
767  C = DRef[i];
768 #endif // ! USE_CL_IN_SHELL
769 
770  B_overline_i[i].MatTMatMul(Ktm, C);
771 // B_overline_i[i].MatTMatMul(Ktm, C);
772 // FullMatrixHandler Gamma(24, 24);
773 // for (integer n = 0; n < NUMNODES; n++) {
774 // Mat3x3 nGamma = Mat3x3(CGR_Rot::MatG, pNode[n]->GetgCurr());
775 // InsertMatrix(Gamma, 1 + 6 * n, 1 + 6 * n, Mat3x3(1.));
776 // InsertMatrix(Gamma, 4 + 6 * n, 4 + 6 * n, nGamma);
777 // InsertMatrix(Gamma, 4 + 6 * n, 4 + 6 * n, Mat3x3(1.));
778 // }
779 
780 
781  Ktm.MatMatMul(Km, B_overline_i[i]);
782 // Ktm.MatMatMul(Kttm, B_overline_i[i]);
783 // Kttm.MatMatMul(Km, Gamma);
784 
785  P_i[i].MatTMatMul(Ktbetaq, C);
786  Ktbetaq.MatMatMul(K_beta_q, B_overline_i[i]);
787 
788  Ktbetaq.MatMatMul(K_beta_beta, P_i[i]);
789 
790 
791 // {
792 // Vec3 n1, n2, m1, m2;
793 // ExtractVec3(n1, stress_i[i], 1);
794 // ExtractVec3(n2, stress_i[i], 4);
795 // ExtractVec3(m1, stress_i[i], 7);
796 // ExtractVec3(m2, stress_i[i], 10);
797 // n1(3) = 0.;
798 // n2(3) = 0.;
799 //
800 // Vec3 Tn1 = T_i[i] * n1;
801 // Vec3 Tn2 = T_i[i] * n2;
802 // Vec3 Tm1 = T_i[i] * m1;
803 // Vec3 Tm2 = T_i[i] * m2;
804 // for (int n = 0; n < NUMNODES; n++) {
805 // for (int m = 0; m < NUMNODES; m++) {
806 // // forze
807 // WM.Add(4 + 6 * n, 4 + 6 * m,
808 // // 1
809 // Phi_Delta_i[i][n].MulTM(Mat3x3(y_i_1[i])) *
810 // Mat3x3(Tn1) * Phi_Delta_i[i][m] *
811 // LI[m](xi_i[i]) * LI[n](xi_i[i])
812 // * alpha_i[i] * w_i[i] * dCoef
813 // +
814 // // 2
815 // Phi_Delta_i[i][n].MulTM(Mat3x3(y_i_2[i])) *
816 // Mat3x3(Tn2) * Phi_Delta_i[i][m] *
817 // LI[m](xi_i[i]) * LI[n](xi_i[i])
818 // * alpha_i[i] * w_i[i] * dCoef
819 // );
820 // WM.Add(4 + 6 * n, 1 + 6 * m,
821 // // 1
822 // Phi_Delta_i[i][n].MulTM(Mat3x3(Tn1)) *
823 // L_alpha_beta_i[i](m + 1, 1) * LI[n](xi_i[i])
824 // * alpha_i[i] * w_i[i] * dCoef
825 // +
826 // // 2
827 // Phi_Delta_i[i][n].MulTM(Mat3x3(Tn2)) *
828 // L_alpha_beta_i[i](m + 1, 2) * LI[n](xi_i[i])
829 // * alpha_i[i] * w_i[i] * dCoef
830 // );
831 // WM.Sub(1 + 6 * n, 4 + 6 * m,
832 // // 1
833 // Mat3x3(Tn1) * Phi_Delta_i[i][m] *
834 // LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 1)
835 // * alpha_i[i] * w_i[i] * dCoef
836 // +
837 // // 2
838 // Mat3x3(Tn2) * Phi_Delta_i[i][m] *
839 // LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 2)
840 // * alpha_i[i] * w_i[i] * dCoef
841 // );
842 // // pezzo in Elle
843 // WM.Add(4 + 6 * n, 4 + 6 * m,
844 // // 1
845 // -Phi_Delta_i[i][n] *
846 // RotManip::Elle(
847 // -phi_tilde_i[i],
848 // T_overline * (Tn1 * Mat3x3(y_i_1[i]))
849 // ) * Phi_Delta_i[i][m]
850 // * LI[n](xi_i[i]) * LI[m](xi_i[i])
851 // * alpha_i[i] * w_i[i] * dCoef
852 // // 2
853 // -
854 // Phi_Delta_i[i][n] *
855 // RotManip::Elle(
856 // -phi_tilde_i[i],
857 // T_overline * (Tn2 * Mat3x3(y_i_2[i]))
858 // ) * Phi_Delta_i[i][m]
859 // * LI[n](xi_i[i]) * LI[m](xi_i[i])
860 // * alpha_i[i] * w_i[i] * dCoef
861 // // 1
862 // +
863 // T_overline.MulMT(RotManip::DRot_IT(phi_tilde_n[n])) *
864 // RotManip::Elle(
865 // -phi_tilde_n[n],
866 // RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
867 // RotManip::DRot(phi_tilde_i[i])
868 // ).MulMT(T_overline) * Tn1 * Mat3x3(y_i_1[i])
869 // ) * Phi_Delta_i[i][m]
870 // * LI[n](xi_i[i]) * LI[m](xi_i[i])
871 // * alpha_i[i] * w_i[i] * dCoef
872 // // 2
873 // +
874 // T_overline.MulMT(RotManip::DRot_IT(phi_tilde_n[n])) *
875 // RotManip::Elle(
876 // -phi_tilde_n[n],
877 // RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
878 // RotManip::DRot(phi_tilde_i[i])
879 // ).MulMT(T_overline) * Tn2 * Mat3x3(y_i_2[i])
880 // ) * Phi_Delta_i[i][m]
881 // * LI[n](xi_i[i]) * LI[m](xi_i[i])
882 // * alpha_i[i] * w_i[i] * dCoef
883 // );
884 //
885 // // momenti
886 // WM.Add(4 + 6 * n, 4 + 6 * m,
887 // // 1
888 // Phi_Delta_i[i][n].MulTM(Mat3x3(k_1_i[i])) * Mat3x3(Tm1) *
889 // Phi_Delta_i[i][m] * LI[m](xi_i[i]) * LI[n](xi_i[i])
890 // * alpha_i[i] * w_i[i] * dCoef
891 // +
892 // Phi_Delta_i[i][n].MulTM(Mat3x3(Tm1)) *
893 // Phi_Delta_i[i][m] * L_alpha_beta_i[i](m + 1, 1)
894 // * LI[n](xi_i[i])
895 // * alpha_i[i] * w_i[i] * dCoef
896 //
897 // +
898 //
899 // // 2
900 // Phi_Delta_i[i][n].MulTM(Mat3x3(k_2_i[i])) * Mat3x3(Tm2) *
901 // Phi_Delta_i[i][m] * LI[m](xi_i[i])
902 // * LI[n](xi_i[i])
903 // * alpha_i[i] * w_i[i] * dCoef
904 // +
905 // Phi_Delta_i[i][n].MulTM(Mat3x3(Tm2)) *
906 // Phi_Delta_i[i][m] * L_alpha_beta_i[i](m + 1, 2)
907 // * LI[n](xi_i[i])
908 // * alpha_i[i] * w_i[i] * dCoef
909 // );
910 // // pezzo in Elle
911 // WM.Add(4 + 6 * n, 4 + 6 * m,
912 // // 1
913 // -
914 // Phi_Delta_i[i][m].MulTMT(
915 // RotManip::Elle(
916 // phi_tilde_i[i],
917 // T_overline.MulTV(Tm1)
918 // )
919 // ) *
920 // RotManip::DRot_I(phi_tilde_n[n]).MulMT(
921 // T_overline
922 // )
923 // * LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 1)
924 // * alpha_i[i] * w_i[i] * dCoef
925 // -
926 // T_overline * RotManip::DRot_IT(phi_tilde_n[m]) *
927 // RotManip::Elle(
928 // phi_tilde_i[i],
929 // T_overline.MulTV(Tm1)
930 // ) * Phi_Delta_i[i][n]
931 // * LI[n](xi_i[i]) * L_alpha_beta_i[i](m + 1, 1)
932 // * alpha_i[i] * w_i[i] * dCoef
933 // // 2
934 // -
935 // Phi_Delta_i[i][m].MulTMT(
936 // RotManip::Elle(
937 // phi_tilde_i[i],
938 // T_overline.MulTV(Tm2)
939 // )
940 // ) *
941 // RotManip::DRot_I(phi_tilde_n[n]).MulMT(
942 // T_overline
943 // )
944 // * LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 2)
945 // * alpha_i[i] * w_i[i] * dCoef
946 // -
947 // T_overline * RotManip::DRot_IT(phi_tilde_n[m]) *
948 // RotManip::Elle(
949 // phi_tilde_i[i],
950 // T_overline.MulTV(Tm2)
951 // ) * Phi_Delta_i[i][n]
952 // * LI[n](xi_i[i]) * L_alpha_beta_i[i](m + 1, 2)
953 // * alpha_i[i] * w_i[i] * dCoef
954 // );
955 // }
956 // WM.Add(4 + 6 * n, 4 + 6 * n,
957 // // 1
958 // T_overline * RotManip::DRot_IT(phi_tilde_n[n]) *
959 // RotManip::Elle(
960 // phi_tilde_n[n],
961 // RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
962 // RotManip::DRot(phi_tilde_i[i])
963 // ).MulMT(T_overline) * Tm1
964 // ) * RotManip::DRot_I(phi_tilde_n[n])
965 // * L_alpha_beta_i[i](n + 1, 1)
966 // * alpha_i[i] * w_i[i] * dCoef
967 // // 2
968 // +
969 // T_overline * RotManip::DRot_IT(phi_tilde_n[n]) *
970 // RotManip::Elle(
971 // phi_tilde_n[n],
972 // RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
973 // RotManip::DRot(phi_tilde_i[i])
974 // ).MulMT(T_overline) * Tm2
975 // ) * RotManip::DRot_I(phi_tilde_n[n])
976 // * L_alpha_beta_i[i](n + 1, 2)
977 // * alpha_i[i] * w_i[i] * dCoef
978 // );
979 // }
980 // }
981 
982 
984 //
985 // VERSIONE DI PIERO
986 //
988  {
989  Vec3 n1, n2, m1, m2;
990  ExtractVec3(n1, stress_i[i], 1);
991  ExtractVec3(n2, stress_i[i], 4);
992  ExtractVec3(m1, stress_i[i], 7);
993  ExtractVec3(m2, stress_i[i], 10);
994 // n1(3) = 0.;
995 // n2(3) = 0.;
996 
997  Vec3 Tn1 = T_i[i] * n1;
998  Vec3 Tn2 = T_i[i] * n2;
999  Vec3 Tm1 = T_i[i] * m1;
1000  Vec3 Tm2 = T_i[i] * m2;
1001  for (int n = 0; n < NUMNODES; n++) {
1002  Mat3x3 ToverGammaITn(T_overline * RotManip::DRot_IT(phi_tilde_n[n]));
1003  for (int m = 0; m < NUMNODES; m++) {
1004  Mat3x3 ToverGammaITm(T_overline * RotManip::DRot_IT(phi_tilde_n[m]));
1005  // forze
1006  WM.Add(4 + 6 * n, 4 + 6 * m,
1007  // 1
1008  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, y_i_1[i])) *
1009  Tn1.Cross(Phi_Delta_i[i][m]) *
1010  LI[m](xi_i[i]) * LI[n](xi_i[i])
1011  * alpha_i[i] * w_i[i] * dCoef
1012  +
1013  // 2
1014  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, y_i_2[i])) *
1015  Tn2.Cross(Phi_Delta_i[i][m]) *
1016  LI[m](xi_i[i]) * LI[n](xi_i[i])
1017  * alpha_i[i] * w_i[i] * dCoef
1018  );
1019  WM.Add(4 + 6 * n, 1 + 6 * m,
1020  // 1
1021  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, Tn1)) *
1022  L_alpha_beta_i[i](m + 1, 1) * LI[n](xi_i[i])
1023  * alpha_i[i] * dCoef
1024  +
1025  // 2
1026  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, Tn2)) *
1027  L_alpha_beta_i[i](m + 1, 2) * LI[n](xi_i[i])
1028  * alpha_i[i] * w_i[i] * w_i[i] * dCoef
1029  );
1030  WM.Sub(1 + 6 * n, 4 + 6 * m,
1031  // 1
1032  Tn1.Cross(Phi_Delta_i[i][m]) *
1033  LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 1)
1034  * alpha_i[i] * w_i[i] * dCoef
1035  +
1036  // 2
1037  Tn2.Cross(Phi_Delta_i[i][m]) *
1038  LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 2)
1039  * alpha_i[i] * w_i[i] * dCoef
1040  );
1041  // pezzo in Elle
1042  WM.Add(4 + 6 * n, 4 + 6 * m,
1043  // 1
1044  -ToverGammaITn *
1045  RotManip::Elle(
1046  -phi_tilde_i[i],
1047  T_overline.MulTM(Mat3x3(MatCross, Tn1)) * y_i_1[i]
1048  ).MulMT(ToverGammaITm)
1049  * LI[n](xi_i[i]) * LI[m](xi_i[i])
1050  * alpha_i[i] * w_i[i] * dCoef
1051  // 2
1052  -ToverGammaITn *
1053  RotManip::Elle(
1054  -phi_tilde_i[i],
1055  T_overline.MulTM(Mat3x3(MatCross, Tn2)) * y_i_2[i]
1056  ).MulMT(ToverGammaITm)
1057  * LI[n](xi_i[i]) * LI[m](xi_i[i])
1058  * alpha_i[i] * w_i[i] * dCoef
1059  );
1060 
1061  // momenti
1062  WM.Add(4 + 6 * n, 4 + 6 * m,
1063  // 1
1064  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, k_1_i[i])) *
1065  Tm1.Cross(Phi_Delta_i[i][m]) * LI[m](xi_i[i]) * LI[n](xi_i[i])
1066  * alpha_i[i] * w_i[i] * dCoef
1067  +
1068  // 2
1069  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, k_2_i[i])) *
1070  Tm2.Cross(Phi_Delta_i[i][m]) * LI[m](xi_i[i])
1071  * LI[n](xi_i[i])
1072  * alpha_i[i] * w_i[i] * dCoef
1073  +
1074  // 1
1075  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, Tm1)) *
1076  Kappa_delta_i_1[i][m]
1077  * LI[n](xi_i[i])
1078  * alpha_i[i] * dCoef
1079  // 2
1080  +
1081  Phi_Delta_i[i][n].MulTM(Mat3x3(MatCross, Tm2)) *
1082  Kappa_delta_i_2[i][m]
1083  * LI[n](xi_i[i])
1084  * alpha_i[i] * w_i[i] * dCoef
1085  );
1086  // pezzo in Elle
1087  WM.Add(4 + 6 * n, 4 + 6 * m,
1088  // 1
1089  -
1090  Phi_Delta_i[i][m].MulTMT(
1091  RotManip::Elle(
1092  phi_tilde_i[i],
1093  T_overline.MulTV(Tm1)
1094  )
1095  ).MulMT(ToverGammaITn)
1096  * LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 1)
1097  * alpha_i[i] * w_i[i] * dCoef
1098  // 2
1099  -
1100  Phi_Delta_i[i][m].MulTMT(
1101  RotManip::Elle(
1102  phi_tilde_i[i],
1103  T_overline.MulTV(Tm2)
1104  )
1105  ).MulMT(ToverGammaITn)
1106  * LI[m](xi_i[i]) * L_alpha_beta_i[i](n + 1, 2)
1107  * alpha_i[i] * w_i[i] * dCoef
1108  // 1
1109  -
1110  ToverGammaITm *
1111  RotManip::Elle(
1112  phi_tilde_i[i],
1113  T_overline.MulTV(Tm1)
1114  )
1115  * Phi_Delta_i[i][n]
1116  * LI[n](xi_i[i]) * L_alpha_beta_i[i](m + 1, 1)
1117  * alpha_i[i] * w_i[i] * dCoef
1118  // 2
1119  -
1120  ToverGammaITm *
1121  RotManip::Elle(
1122  phi_tilde_i[i],
1123  T_overline.MulTV(Tm2)
1124  )
1125  * Phi_Delta_i[i][n]
1126  * LI[n](xi_i[i]) * L_alpha_beta_i[i](m + 1, 2)
1127  * alpha_i[i] * w_i[i] * dCoef
1128  );
1129  }
1130  //pezzo in Elle sono su n per le forze
1131  WM.Add(4 + 6 * n, 4 + 6 * n,
1132  // 1
1133  ToverGammaITn *
1134  RotManip::Elle(
1135  -phi_tilde_n[n],
1136  RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
1137  RotManip::DRot(phi_tilde_i[i])
1138  ).MulMT(T_overline) * Tn1.Cross(y_i_1[i])
1139  ).MulMT(ToverGammaITn)
1140  * LI[n](xi_i[i])
1141  * alpha_i[i] * w_i[i] * dCoef
1142  // 2
1143  +
1144  ToverGammaITn *
1145  RotManip::Elle(
1146  -phi_tilde_n[n],
1147  RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
1148  RotManip::DRot(phi_tilde_i[i])
1149  ).MulMT(T_overline) * Tn2.Cross(y_i_2[i])
1150  ).MulMT(ToverGammaITn)
1151  * LI[n](xi_i[i])
1152  * alpha_i[i] * w_i[i] * dCoef
1153  );
1154  //pezzo in Elle sono su n per i momenti
1155  WM.Add(4 + 6 * n, 4 + 6 * n,
1156  // 1
1157  ToverGammaITn *
1158  RotManip::Elle(
1159  phi_tilde_n[n],
1160  RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
1161  RotManip::DRot(phi_tilde_i[i])
1162  ).MulMT(T_overline) * Tm1
1163  ).MulMT(ToverGammaITn)
1164  * L_alpha_beta_i[i](n + 1, 1)
1165  * alpha_i[i] * w_i[i] * dCoef
1166  // 2
1167  +
1168  ToverGammaITn *
1169  RotManip::Elle(
1170  phi_tilde_n[n],
1171  RotManip::DRot_IT(phi_tilde_n[n]).MulMT(
1172  RotManip::DRot(phi_tilde_i[i])
1173  ).MulMT(T_overline) * Tm2
1174  ).MulMT(ToverGammaITn)
1175  * L_alpha_beta_i[i](n + 1, 2)
1176  * alpha_i[i] * w_i[i] * dCoef
1177  );
1178  }
1179  }
1180 
1181 
1182 
1183 // std::cerr << "Kg:\n" << std::fixed << std::setprecision(12) << Kg << std::endl;
1184 
1185 
1186 // std::cerr << "Km:\n" << std::fixed << std::setprecision(12) << Km << std::endl;
1187  // AssembleMatrix(WM, 1, 1, Km, alpha_i[i] * w_i[i] * dCoef);
1188 
1189  // AssembleTransposeMatrix(WM, 1, 25, K_beta_q, alpha_i[i] * w_i[i]);
1190  // AssembleMatrix(WM, 25, 1, K_beta_q, alpha_i[i] * w_i[i] * dCoef);
1191  // AssembleMatrix(WM, 25, 25, K_beta_beta, alpha_i[i] * w_i[i]);
1192 
1193  WM.Add(1, 1, Km, alpha_i[i] * w_i[i] * dCoef);
1194 
1195  WM.AddT(1, 25, K_beta_q, alpha_i[i] * w_i[i]);
1196  WM.Add(25, 1, K_beta_q, alpha_i[i] * w_i[i] * dCoef);
1197  WM.Add(25, 25, K_beta_beta, alpha_i[i] * w_i[i]);
1198  }
1199  return WorkMat;
1200 
1201 }
1202 
1203 // Contribution to restart file
1204 std::ostream&
1205 Shell4EAS::Restart(std::ostream& out) const
1206 {
1207  return out << "# not implemented yet" << std::endl;
1208 }
1209 
1210 // Initial settings
1211 void
1212 Shell4EAS::SetValue(DataManager *pDM,
1213  VectorHandler& /* X */ , VectorHandler& /* XP */ ,
1214  SimulationEntity::Hints *ph)
1215 {
1216  NO_OP;
1217 }
1218 
1219 VariableSubMatrixHandler&
1220 Shell4EAS::InitialAssJac(VariableSubMatrixHandler& WorkMat,
1221  const VectorHandler& XCurr)
1222 {
1223  WorkMat.SetNullMatrix();
1224  return WorkMat;
1225 }
1226 
1227 SubVectorHandler&
1228 Shell4EAS::InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr)
1229 {
1230  WorkVec.Resize(0);
1231  return WorkVec;
1232 }
1233 
1234 // Access to nodes
1235 const StructNode*
1236 Shell4EAS::pGetNode(unsigned int i) const
1237 {
1238  ASSERT(i >= 1 && i <= 4);
1239  switch (i) {
1240  case 1:
1241  case 2:
1242  case 3:
1243  case 4:
1244  return pNode[i - 1];
1245  default:
1246  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1247  }
1248 }
1249 
1250 void
1251 Shell4EAS::Output(OutputHandler& OH) const
1252 {
1253  if (bToBeOutput()) {
1254  if (OH.UseText(OutputHandler::PLATES)) {
1255  std::ostream& out = OH.Plates();
1256  out << std::setw(8) << GetLabel();
1257  // TODO: complete
1258  for (integer i = 0; i < NUMIP; i++) {
1259  for (integer r = 1; r <= 12; r++) {
1260  out << " " << stress_i[i](r);
1261  }
1262  }
1263  out << std::endl;
1264  }
1265  }
1266 }
1267 
1268 
1269 typedef LinearElasticGenericConstitutiveLaw<Shell::vh, Shell::fmh> LEGCLShell;
1270 
1271 Elem*
1272 ReadShell4EAS(DataManager* pDM,
1273  MBDynParser& HP,
1274  const DofOwner* pDO,
1275  unsigned int uLabel)
1276 {
1277  silent_cout("Shell4EAS(" << uLabel << "): warning, \"shell4eas\" is deprecated; use \"shell4easans\" instead" << std::endl);
1278 
1279  const StructNode* pN[4];
1280  Mat3x3 R[4];
1281  for (unsigned i = 0; i < 4; i++) {
1282  pN[i] = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1283  ReferenceFrame RF(pN[i]);
1284  if (HP.IsKeyWord("orientation")) {
1285  R[i] = HP.GetRotRel(RF);
1286  } else {
1287  R[i] = Eye3;
1288  }
1289  }
1290 
1291  /* Prestress and prestrain */
1292  Shell::vh PreStress(12);
1293  PreStress.Reset();
1294 
1295 #ifdef USE_CL_IN_SHELL
1296  const ConstitutiveLaw<Shell::vh, Shell::fmh> *pD[4];
1297 
1298  TplDriveCaller<Shell::vh>* pTplDC = new ZeroTplDriveCaller<Shell::vh>;
1299 
1300  Shell::fmh S(12, 12);
1301  for (unsigned ir = 1; ir <= 12; ir++) {
1302  for (unsigned ic = 1; ic <= 12; ic++) {
1303  S(ir, ic) = HP.GetReal();
1304  }
1305  }
1306 
1307  pD[0] = 0;
1308  SAFENEWWITHCONSTRUCTOR(pD[0], LEGCLShell, LEGCLShell(pTplDC, PreStress, S));
1309 
1310  for (unsigned i = 1; i < 4; i++) {
1311  pD[i] = pD[0]->Copy();
1312  }
1313 #else // ! USE_CL_IN_SHELL
1314  Shell::fmh pD(12, 12);
1315 
1316  if (ReadShellConstLaw(HP, pD, PreStress)) {
1317  silent_cerr("Shell(" << uLabel << "): unable to read constitutive law" << std::endl);
1318  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1319  }
1320 #endif // ! USE_CL_IN_SHELL
1321 
1322  flag fOut = pDM->fReadOutput(HP, Elem::PLATE);
1323 
1324  Elem *pEl = 0;
1325  SAFENEWWITHCONSTRUCTOR(pEl, Shell4EAS,
1326  Shell4EAS(uLabel, pDO,
1327  pN[0], pN[1], pN[2], pN[3],
1328  R[0], R[1], R[2], R[3],
1329 #ifdef USE_CL_IN_SHELL
1330 #else // ! USE_CL_IN_SHELL
1331  pD, PreStress,
1332 #endif // ! USE_CL_IN_SHELL
1333  fOut));
1334 
1335  std::ostream& out = pDM->GetLogFile();
1336  out << "shell4: " << uLabel
1337  << " " << pN[0]->GetLabel()
1338  << " " << pN[1]->GetLabel()
1339  << " " << pN[2]->GetLabel()
1340  << " " << pN[3]->GetLabel()
1341  << std::endl;
1342 
1343  return pEl;
1344 }
1345 
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
constltp_impl.h
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
Shell4EAS::T_overline
Mat3x3 T_overline
Definition: shelleas.h:161
Shell4EAS::ComputeInitialNodeAndIptOrientation
void ComputeInitialNodeAndIptOrientation()
Definition: shelleas.cc:123
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Shell4EAS::Shell4EAS
Shell4EAS(unsigned int uL, const DofOwner *pDO, const StructNode *pN1, const StructNode *pN2, const StructNode *pN3, const StructNode *pN4, const Mat3x3 &R1, const Mat3x3 &R2, const Mat3x3 &R3, const Mat3x3 &R4, const fmh &pDTmp, const vh &PreStress, flag fOut)
Definition: shelleas.cc:215
Shell4EAS::P_i
vfmh P_i
Definition: shelleas.h:230
RotManip::DRot_IT
Mat3x3 DRot_IT(const Vec3 &phi)
Definition: Rot.cc:100
Shell4EAS::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: shelleas.cc:1205
Shell4EAS::SetValue
void SetValue(DataManager *pDM, VectorHandler &, VectorHandler &, SimulationEntity::Hints *ph=0)
Definition: shelleas.cc:1212
Zero3
const Vec3 Zero3(0., 0., 0.)
Vec3::Cross
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
Shell4EAS::S_alpha_beta_0
fmh S_alpha_beta_0
Definition: shelleas.h:223
Shell4EAS::epsilon_hat
vh epsilon_hat
Definition: shelleas.h:236
Shell4EAS::AssJac
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: shelleas.cc:545
flag
long int flag
Definition: mbdyn.h:43
StructNode::GetRRef
virtual const Mat3x3 & GetRRef(void) const
Definition: strnode.h:1006
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
Shell4EAS::y_i_1
Vec3 y_i_1[NUMIP]
Definition: shelleas.h:232
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
Shell4EAS::T_0_i
Mat3x3 T_0_i[NUMIP]
Definition: shelleas.h:166
Shell4EAS::k_tilde_1_i
Vec3 k_tilde_1_i[NUMIP]
Definition: shelleas.h:202
Shell4EAS::xi_i
static doublereal xi_i[NUMIP][2]
Definition: shelleas.h:105
VectorHandler::ResizeReset
virtual void ResizeReset(integer)
Definition: vh.cc:55
tpldrive_impl.h
MatCross
const MatCross_Manip MatCross
Definition: matvec3.cc:639
Shell4EAS::y_i_2
Vec3 y_i_2[NUMIP]
Definition: shelleas.h:233
Mat3x3::MulTMT
Mat3x3 MulTMT(const Mat3x3 &m) const
Definition: matvec3.cc:538
Shell4EAS::bPreStress
bool bPreStress
Definition: shelleas.h:244
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
Shell4EAS::L_alpha_beta_i
vfmh L_alpha_beta_i
Definition: shelleas.h:226
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
Shell4EAS::Output
virtual void Output(OutputHandler &OH) const
Definition: shelleas.cc:1251
Vec3::Norm
doublereal Norm(void) const
Definition: matvec3.h:263
Shell4EAS::Phi_Delta_i
Mat3x3 Phi_Delta_i[NUMIP][NUMNODES]
Definition: shelleas.h:171
FullSubMatrixHandler::Add
void Add(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:209
Shell4EAS::stress_i
vvh stress_i
Definition: shelleas.h:252
RotManip::DRot
Mat3x3 DRot(const Vec3 &phi)
Definition: Rot.cc:74
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
Mat3x3::GetCol
Vec3 GetCol(unsigned short int i) const
Definition: matvec3.h:903
Shell4EAS::Kappa_delta_i_1
Mat3x3 Kappa_delta_i_1[NUMIP][NUMNODES]
Definition: shelleas.h:172
Shell4EAS::k_tilde_2_i
Vec3 k_tilde_2_i[NUMIP]
Definition: shelleas.h:203
Shell4EAS::T0_overline
Mat3x3 T0_overline
Definition: shelleas.h:159
Shell4EAS::phi_tilde_i
Vec3 phi_tilde_i[NUMIP]
Definition: shelleas.h:156
MyVectorHandler::Put
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:699
NO_OP
#define NO_OP
Definition: myassert.h:74
Shell4EAS::ComputeIPCurvature
void ComputeIPCurvature()
Definition: shelleas.cc:191
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
Shell4EAS::xi_0
static doublereal xi_0[2]
Definition: shelleas.h:127
Shell4EAS::eps_tilde_1_0_i
Vec3 eps_tilde_1_0_i[NUMIP]
Definition: shelleas.h:191
Shell4EAS::xi_n
static doublereal xi_n[NUMNODES][2]
Definition: shelleas.h:125
Shell4EAS::~Shell4EAS
virtual ~Shell4EAS(void)
Definition: shelleas.cc:410
Shell4EAS::alpha_i
doublereal alpha_i[NUMIP]
Definition: shelleas.h:225
FullSubMatrixHandler::AddT
void AddT(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:227
RotManip::VecRot
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
Shell4EAS::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: shelleas.h:327
SubVectorHandler::PutRowIndex
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
Shell4EAS::T_i
Mat3x3 T_i[NUMIP]
Definition: shelleas.h:169
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
Shell4EAS::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: shelleas.h:297
ReadShellConstLaw
int ReadShellConstLaw(MBDynParser &HP, Shell::fmh &pD, Shell::vh &PreStress)
Definition: shell.cc:62
Shell4EAS::DRef
vfmh DRef
Definition: shelleas.h:249
Shell4EAS::phi_tilde_n
Vec3 phi_tilde_n[NUMNODES]
Definition: shelleas.h:153
Shell4EAS::k_tilde_1_0_i
Vec3 k_tilde_1_0_i[NUMIP]
Definition: shelleas.h:199
Shell4EAS::InitialAssRes
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: shelleas.cc:1228
Shell4EAS::Kappa_delta_i_2
Mat3x3 Kappa_delta_i_2[NUMIP][NUMNODES]
Definition: shelleas.h:173
VariableSubMatrixHandler::SetNullMatrix
void SetNullMatrix(void)
Definition: submat.h:1159
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
Shell4EAS::xa
Vec3 xa[NUMNODES]
Definition: shelleas.h:148
Shell4EAS::AssRes
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: shelleas.cc:422
Shell4EAS::w_i
static doublereal w_i[NUMIP]
Definition: shelleas.h:106
RotManip::Rot
Mat3x3 Rot(const Vec3 &phi)
Definition: Rot.cc:62
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
Shell4EAS::alpha_0
doublereal alpha_0
Definition: shelleas.h:224
StructDispNode::iGetFirstPositionIndex
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
Shell4EAS::k_1_i
Vec3 k_1_i[NUMIP]
Definition: shelleas.h:184
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
ReadShell4EAS
Elem * ReadShell4EAS(DataManager *pDM, MBDynParser &HP, const DofOwner *pDO, unsigned int uLabel)
Definition: shelleas.cc:1272
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
Mat3x3::MulTM
Mat3x3 MulTM(const Mat3x3 &m) const
Definition: matvec3.cc:500
Shell4EAS::epsilon
vh epsilon
Definition: shelleas.h:237
MyVectorHandler::Reset
virtual void Reset(void)
Definition: vh.cc:459
Shell4EAS::PreStress
vh PreStress
Definition: shelleas.h:245
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
RotManip::Elle
Mat3x3 Elle(const Vec3 &phi, const Vec3 &a)
Definition: Rot.cc:179
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
Shell4EAS::iTa_i
Mat3x3 iTa_i[NUMIP]
Definition: shelleas.h:151
grad::Norm
VectorExpression< VectorExpr, N_rows >::ScalarType Norm(const VectorExpression< VectorExpr, N_rows > &u)
Definition: matvec.h:3163
Shell4EAS::iTa
Mat3x3 iTa[NUMNODES]
Definition: shelleas.h:150
Shell4EAS::eps_tilde_2_0_i
Vec3 eps_tilde_2_0_i[NUMIP]
Definition: shelleas.h:192
Shell4EAS::B_overline_i
vfmh B_overline_i
Definition: shelleas.h:228
Elem
Definition: elem.h:75
Shell4EAS::xa_0
Vec3 xa_0[NUMNODES]
Definition: shelleas.h:147
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
RotManip::RotAndDRot
void RotAndDRot(const Vec3 &phi, Mat3x3 &Phi, Mat3x3 &Ga)
Definition: Rot.cc:86
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
FullSubMatrixHandler::Sub
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
a
static const doublereal a
Definition: hfluid_.h:289
Shell4EAS::pGetNode
virtual const StructNode * pGetNode(unsigned int i) const
Definition: shelleas.cc:1236
Shell4EAS::UpdateNodalAndAveragePosAndOrientation
void UpdateNodalAndAveragePosAndOrientation()
Definition: shelleas.cc:101
Mat3x3::MulMT
Mat3x3 MulMT(const Mat3x3 &m) const
Definition: matvec3.cc:444
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
Shell4EAS::T_0_0
Mat3x3 T_0_0
Definition: shelleas.h:165
OutputHandler::Plates
std::ostream & Plates(void) const
Definition: output.h:580
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
LEGCLShell
LinearElasticGenericConstitutiveLaw< Shell::vh, Shell::fmh > LEGCLShell
Definition: shelleas.cc:1269
MatrixHandler::MatMatMul
virtual MatrixHandler & MatMatMul(MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:232
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
MatrixHandler::MatTMatMul
virtual MatrixHandler & MatTMatMul(MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:241
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
Shell4EAS::eps_tilde_1_i
Vec3 eps_tilde_1_i[NUMIP]
Definition: shelleas.h:194
Shell4EAS::pNode
const StructNode * pNode[NUMNODES]
Definition: shelleas.h:138
Shell4EAS::T_0
Mat3x3 T_0
Definition: shelleas.h:168
Shell4EAS::k_2_i
Vec3 k_2_i[NUMIP]
Definition: shelleas.h:185
Shell::ConstitutiveLawOwnerType
ConstitutiveLawOwner< vh, fmh > ConstitutiveLawOwnerType
Definition: shell.h:84
Shell4EAS::eps_tilde_2_i
Vec3 eps_tilde_2_i[NUMIP]
Definition: shelleas.h:195
Shell4EAS::InitialAssJac
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: shelleas.cc:1220
Shell4EAS::InterpolateOrientation
void InterpolateOrientation()
Definition: shelleas.cc:157
VectorHandler::Resize
virtual void Resize(integer iNewSize)=0
RotManip::DRot_I
Mat3x3 DRot_I(const Vec3 &phi)
Definition: Rot.cc:111
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
OutputHandler::UseText
bool UseText(int out) const
Definition: output.cc:446
Shell
Definition: shell.h:63
Shell4EAS::k_tilde_2_0_i
Vec3 k_tilde_2_0_i[NUMIP]
Definition: shelleas.h:200
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056