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: Bryan Peterson, Antonio Recuero, Radu Serban
13 // =============================================================================
14 // Hexahedronal element with 8 nodes (with EAS)
15 // =============================================================================
16
17 //// RADU
18 //// A lot more to do here...
19 //// - reconsider the use of large static matrices
20 //// - more use of Eigen expressions
21
22 #include "chrono/core/ChException.h"
23 #include "chrono/physics/ChSystem.h"
24 #include "chrono/fea/ChElementHexaANCF_3813.h"
25
26 namespace chrono {
27 namespace fea {
28
29 // -----------------------------------------------------------------------------
30
ChElementHexaANCF_3813()31 ChElementHexaANCF_3813::ChElementHexaANCF_3813() : m_flag_HE(ANALYTICAL) {
32 m_nodes.resize(8);
33 }
34
35 // -----------------------------------------------------------------------------
SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA,std::shared_ptr<ChNodeFEAxyz> nodeB,std::shared_ptr<ChNodeFEAxyz> nodeC,std::shared_ptr<ChNodeFEAxyz> nodeD,std::shared_ptr<ChNodeFEAxyz> nodeE,std::shared_ptr<ChNodeFEAxyz> nodeF,std::shared_ptr<ChNodeFEAxyz> nodeG,std::shared_ptr<ChNodeFEAxyz> nodeH)36 void ChElementHexaANCF_3813::SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA,
37 std::shared_ptr<ChNodeFEAxyz> nodeB,
38 std::shared_ptr<ChNodeFEAxyz> nodeC,
39 std::shared_ptr<ChNodeFEAxyz> nodeD,
40 std::shared_ptr<ChNodeFEAxyz> nodeE,
41 std::shared_ptr<ChNodeFEAxyz> nodeF,
42 std::shared_ptr<ChNodeFEAxyz> nodeG,
43 std::shared_ptr<ChNodeFEAxyz> nodeH) {
44 assert(nodeA);
45 assert(nodeB);
46 assert(nodeC);
47 assert(nodeD);
48 assert(nodeE);
49 assert(nodeF);
50 assert(nodeG);
51 assert(nodeH);
52
53 m_nodes[0] = nodeA;
54 m_nodes[1] = nodeB;
55 m_nodes[2] = nodeC;
56 m_nodes[3] = nodeD;
57 m_nodes[4] = nodeE;
58 m_nodes[5] = nodeF;
59 m_nodes[6] = nodeG;
60 m_nodes[7] = nodeH;
61 std::vector<ChVariables*> mvars;
62 mvars.push_back(&m_nodes[0]->Variables());
63 mvars.push_back(&m_nodes[1]->Variables());
64 mvars.push_back(&m_nodes[2]->Variables());
65 mvars.push_back(&m_nodes[3]->Variables());
66 mvars.push_back(&m_nodes[4]->Variables());
67 mvars.push_back(&m_nodes[5]->Variables());
68 mvars.push_back(&m_nodes[6]->Variables());
69 mvars.push_back(&m_nodes[7]->Variables());
70 Kmatr.SetVariables(mvars);
71 // EAS
72 // Initial position
73 ChVector<> pA = m_nodes[0]->GetPos();
74 ChVector<> pB = m_nodes[1]->GetPos();
75 ChVector<> pC = m_nodes[2]->GetPos();
76 ChVector<> pD = m_nodes[3]->GetPos();
77 ChVector<> pE = m_nodes[4]->GetPos();
78 ChVector<> pF = m_nodes[5]->GetPos();
79 ChVector<> pG = m_nodes[6]->GetPos();
80 ChVector<> pH = m_nodes[7]->GetPos();
81 m_d0(0, 0) = pA[0];
82 m_d0(0, 1) = pA[1];
83 m_d0(0, 2) = pA[2];
84 m_d0(1, 0) = pB[0];
85 m_d0(1, 1) = pB[1];
86 m_d0(1, 2) = pB[2];
87 m_d0(2, 0) = pC[0];
88 m_d0(2, 1) = pC[1];
89 m_d0(2, 2) = pC[2];
90 m_d0(3, 0) = pD[0];
91 m_d0(3, 1) = pD[1];
92 m_d0(3, 2) = pD[2];
93 m_d0(4, 0) = pE[0];
94 m_d0(4, 1) = pE[1];
95 m_d0(4, 2) = pE[2];
96 m_d0(5, 0) = pF[0];
97 m_d0(5, 1) = pF[1];
98 m_d0(5, 2) = pF[2];
99 m_d0(6, 0) = pG[0];
100 m_d0(6, 1) = pG[1];
101 m_d0(6, 2) = pG[2];
102 m_d0(7, 0) = pH[0];
103 m_d0(7, 1) = pH[1];
104 m_d0(7, 2) = pH[2];
105 // EAS
106 }
107
108 // -----------------------------------------------------------------------------
109
SetStockAlpha(double a1,double a2,double a3,double a4,double a5,double a6,double a7,double a8,double a9)110 void ChElementHexaANCF_3813::SetStockAlpha(double a1,
111 double a2,
112 double a3,
113 double a4,
114 double a5,
115 double a6,
116 double a7,
117 double a8,
118 double a9) {
119 m_stock_alpha_EAS(0) = a1;
120 m_stock_alpha_EAS(1) = a2;
121 m_stock_alpha_EAS(2) = a3;
122 m_stock_alpha_EAS(3) = a4;
123 m_stock_alpha_EAS(4) = a5;
124 m_stock_alpha_EAS(5) = a6;
125 m_stock_alpha_EAS(6) = a7;
126 m_stock_alpha_EAS(7) = a8;
127 m_stock_alpha_EAS(8) = a9;
128 }
129
130 // -----------------------------------------------------------------------------
131
132 // Internal force, EAS stiffness, and analytical jacobian are calculated
133 class Brick_ForceAnalytical : public ChIntegrable3D<ChVectorN<double, 906>> {
134 public:
135 Brick_ForceAnalytical(ChMatrixNM<double, 8, 3>* d_,
136 ChMatrixNM<double, 8, 3>* d0_,
137 ChElementHexaANCF_3813* element_,
138 ChMatrixNM<double, 6, 6>* T0_,
139 double* detJ0C_,
140 ChVectorN<double, 9>* alpha_eas_);
141
142 Brick_ForceAnalytical(ChMatrixNM<double, 8, 3>* d_,
143 ChMatrixNM<double, 8, 3>* d0_,
144 ChElementHexaANCF_3813* element_,
145 ChMatrixNM<double, 6, 6>* T0_,
146 double* detJ0C_,
147 ChVectorN<double, 9>* alpha_eas_,
148 double* E_,
149 double* v_);
~Brick_ForceAnalytical()150 ~Brick_ForceAnalytical() {}
151
152 private:
153 ChElementHexaANCF_3813* element;
154 ChMatrixNM<double, 8, 3>* d; // Pointer to a matrix containing the element coordinates
155 ChMatrixNM<double, 8, 3>* d0; // Pointer to a matrix containing the element initial coordinates
156 ChMatrixNM<double, 6, 6>* T0; // Pointer to transformation matrix for Enhanced Assumed Strain (EAS)
157 ChVectorN<double, 9>* alpha_eas; // Pointer to the 9 internal parameters for EAS
158 double* detJ0C; // Pointer to determinant of the initial Jacobian at the element center
159 double* E; // Pointer to Young modulus
160 double* v; // Pointer to Poisson ratio
161
162 ChVectorN<double, 24> Fint; // Generalized internal (elastic) force vector
163 ChMatrixNM<double, 24, 24> JAC11; // Jacobian of internal forces for implicit numerical integration
164 ChMatrixNM<double, 9, 24> Gd; // Jacobian (w.r.t. coordinates) of the initial pos. vector gradient matrix
165 ChVectorN<double, 6> stress; // stress tensor in vector form
166 ChMatrixNM<double, 9, 9> Sigm; // stress tensor in sparse form
167 ChMatrixNM<double, 6, 6> E_eps; // Matrix of elastic coefficients (features orthotropy)
168 ChMatrixNM<double, 3, 24> Sx; // Sparse shape function matrix, X derivative
169 ChMatrixNM<double, 3, 24> Sy; // Sparse shape function matrix, Y derivative
170 ChMatrixNM<double, 3, 24> Sz; // Sparse shape function matrix, Z derivative
171 ChElementHexaANCF_3813::ShapeVector Nx; // Dense shape function vector, X derivative
172 ChElementHexaANCF_3813::ShapeVector Ny; // Dense shape function vector, Y derivative
173 ChElementHexaANCF_3813::ShapeVector Nz; // Dense shape function vector, Z derivative
174 ChMatrixNM<double, 6, 24> strainD; // Derivative of the strains w.r.t. the coordinates. Includes orthotropy
175 ChVectorN<double, 6> strain; // Vector of strains
176 double detJ0; // Determinant of the initial position vector gradient matrix
177 // EAS
178 ChMatrixNM<double, 6, 9> M; // Shape function matrix for Enhanced Assumed Strain
179 ChMatrixNM<double, 6, 9> G; // Matrix G interpolates the internal parameters of EAS
180 ChVectorN<double, 6> strain_EAS; // Enhanced assumed strain vector
181
182 // Evaluate (strainD'*strain) at a point
183 virtual void Evaluate(ChVectorN<double, 906>& result, const double x, const double y, const double z) override;
184 };
185
Brick_ForceAnalytical(ChMatrixNM<double,8,3> * d_,ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_,ChMatrixNM<double,6,6> * T0_,double * detJ0C_,ChVectorN<double,9> * alpha_eas_)186 Brick_ForceAnalytical::Brick_ForceAnalytical(ChMatrixNM<double, 8, 3>* d_,
187 ChMatrixNM<double, 8, 3>* d0_,
188 ChElementHexaANCF_3813* element_,
189 ChMatrixNM<double, 6, 6>* T0_,
190 double* detJ0C_,
191 ChVectorN<double, 9>* alpha_eas_)
192 : element(element_), d(d_), d0(d0_), T0(T0_), alpha_eas(alpha_eas_), detJ0C(detJ0C_) {
193 E_eps.setZero();
194 Gd.setZero();
195 Sigm.setZero();
196
197 Sx.setZero();
198 Sy.setZero();
199 Sz.setZero();
200 }
201
Brick_ForceAnalytical(ChMatrixNM<double,8,3> * d_,ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_,ChMatrixNM<double,6,6> * T0_,double * detJ0C_,ChVectorN<double,9> * alpha_eas_,double * E_,double * v_)202 Brick_ForceAnalytical::Brick_ForceAnalytical(ChMatrixNM<double, 8, 3>* d_,
203 ChMatrixNM<double, 8, 3>* d0_,
204 ChElementHexaANCF_3813* element_,
205 ChMatrixNM<double, 6, 6>* T0_,
206 double* detJ0C_,
207 ChVectorN<double, 9>* alpha_eas_,
208 double* E_,
209 double* v_)
210 : element(element_), d(d_), d0(d0_), T0(T0_), alpha_eas(alpha_eas_), detJ0C(detJ0C_), E(E_), v(v_) {
211 E_eps.setZero();
212 Gd.setZero();
213 Sigm.setZero();
214
215 Sx.setZero();
216 Sy.setZero();
217 Sz.setZero();
218 }
219
Evaluate(ChVectorN<double,906> & result,const double x,const double y,const double z)220 void Brick_ForceAnalytical::Evaluate(ChVectorN<double, 906>& result, const double x, const double y, const double z) {
221 element->ShapeFunctionsDerivativeX(Nx, x, y, z);
222 element->ShapeFunctionsDerivativeY(Ny, x, y, z);
223 element->ShapeFunctionsDerivativeZ(Nz, x, y, z);
224
225 element->Basis_M(M, x, y, z); // EAS
226
227 if (!element->m_isMooney) { // m_isMooney == false means use linear material
228 double DD = (*E) * (1.0 - (*v)) / ((1.0 + (*v)) * (1.0 - 2.0 * (*v)));
229 E_eps.fillDiagonal(1.0);
230 E_eps(0, 1) = (*v) / (1.0 - (*v));
231 E_eps(0, 3) = (*v) / (1.0 - (*v));
232 E_eps(1, 0) = (*v) / (1.0 - (*v));
233 E_eps(1, 3) = (*v) / (1.0 - (*v));
234 E_eps(2, 2) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
235 E_eps(3, 0) = (*v) / (1.0 - (*v));
236 E_eps(3, 1) = (*v) / (1.0 - (*v));
237 E_eps(4, 4) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
238 E_eps(5, 5) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
239 E_eps *= DD;
240 }
241
242 // Sd=[Nd1*eye(3) Nd2*eye(3) Nd3*eye(3) Nd4*eye(3)]
243
244 for (int i = 0; i < 8; i++) {
245 Sx(0, 3 * i + 0) = Nx(i);
246 Sx(1, 3 * i + 1) = Nx(i);
247 Sx(2, 3 * i + 2) = Nx(i);
248 }
249
250 for (int i = 0; i < 8; i++) {
251 Sy(0, 3 * i + 0) = Ny(i);
252 Sy(1, 3 * i + 1) = Ny(i);
253 Sy(2, 3 * i + 2) = Ny(i);
254 }
255
256 for (int i = 0; i < 8; i++) {
257 Sz(0, 3 * i + 0) = Nz(i);
258 Sz(1, 3 * i + 1) = Nz(i);
259 Sz(2, 3 * i + 2) = Nz(i);
260 }
261
262 // EAS and Initial Shape
263 ChMatrixNM<double, 3, 3> rd0;
264 rd0.col(0) = (*d0).transpose() * Nx.transpose();
265 rd0.col(1) = (*d0).transpose() * Ny.transpose();
266 rd0.col(2) = (*d0).transpose() * Nz.transpose();
267 detJ0 = rd0.determinant();
268
269 // Transformation : Orthogonal transformation (A and J)
270
271 ChVector<double> G1;
272 ChVector<double> G2;
273 ChVector<double> G3;
274 ChVector<double> G1xG2;
275 G1[0] = rd0(0, 0);
276 G2[0] = rd0(0, 1);
277 G3[0] = rd0(0, 2);
278 G1[1] = rd0(1, 0);
279 G2[1] = rd0(1, 1);
280 G3[1] = rd0(1, 2);
281 G1[2] = rd0(2, 0);
282 G2[2] = rd0(2, 1);
283 G3[2] = rd0(2, 2);
284 G1xG2.Cross(G1, G2);
285
286 // Tangent Frame
287 ChVector<> A1 = G1 / sqrt(G1[0] * G1[0] + G1[1] * G1[1] + G1[2] * G1[2]);
288 ChVector<> A3 = G1xG2 / sqrt(G1xG2[0] * G1xG2[0] + G1xG2[1] * G1xG2[1] + G1xG2[2] * G1xG2[2]);
289 ChVector<> A2 = A3.Cross(A1);
290
291 // Direction for orthotropic material
292 double theta = 0.0;
293 ChVector<> AA1 = A1 * cos(theta) + A2 * sin(theta);
294 ChVector<> AA2 = -A1 * sin(theta) + A2 * cos(theta);
295 ChVector<> AA3 = A3;
296
297 // Beta
298 ChMatrixNM<double, 3, 3> j0 = rd0.inverse();
299 ChVector<double> j01;
300 ChVector<double> j02;
301 ChVector<double> j03;
302 j01[0] = j0(0, 0);
303 j02[0] = j0(1, 0);
304 j03[0] = j0(2, 0);
305 j01[1] = j0(0, 1);
306 j02[1] = j0(1, 1);
307 j03[1] = j0(2, 1);
308 j01[2] = j0(0, 2);
309 j02[2] = j0(1, 2);
310 j03[2] = j0(2, 2);
311 ChVectorN<double, 9> beta;
312 beta(0) = Vdot(AA1, j01);
313 beta(1) = Vdot(AA2, j01);
314 beta(2) = Vdot(AA3, j01);
315 beta(3) = Vdot(AA1, j02);
316 beta(4) = Vdot(AA2, j02);
317 beta(5) = Vdot(AA3, j02);
318 beta(6) = Vdot(AA1, j03);
319 beta(7) = Vdot(AA2, j03);
320 beta(8) = Vdot(AA3, j03);
321
322 // Enhanced Assumed Strain
323 G = (*T0) * M * ((*detJ0C) / (detJ0));
324 strain_EAS = G * (*alpha_eas);
325
326 ChMatrixNM<double, 8, 8> d_d = (*d) * (*d).transpose();
327 ChVectorN<double, 8> ddNx = d_d * Nx.transpose();
328 ChVectorN<double, 8> ddNy = d_d * Ny.transpose();
329 ChVectorN<double, 8> ddNz = d_d * Nz.transpose();
330
331 ChMatrixNM<double, 8, 8> d0_d0 = (*d0) * (*d0).transpose();
332 ChVectorN<double, 8> d0d0Nx = d0_d0 * Nx.transpose();
333 ChVectorN<double, 8> d0d0Ny = d0_d0 * Ny.transpose();
334 ChVectorN<double, 8> d0d0Nz = d0_d0 * Nz.transpose();
335
336 // Strain component
337
338 ChVectorN<double, 6> strain_til;
339 strain_til(0) = 0.5 * Nx.dot(ddNx - d0d0Nx);
340 strain_til(1) = 0.5 * Ny.dot(ddNy - d0d0Ny);
341 strain_til(2) = Nx.dot(ddNy - d0d0Ny);
342 //== Compatible strain (No ANS) ==//
343 strain_til(3) = 0.5 * Nz.dot(ddNz - d0d0Nz);
344 strain_til(4) = Nx.dot(ddNz - d0d0Nz);
345 strain_til(5) = Ny.dot(ddNz - d0d0Nz);
346
347 // //// For orthotropic material ///
348 strain(0) = strain_til(0) * beta(0) * beta(0) + strain_til(1) * beta(3) * beta(3) +
349 strain_til(2) * beta(0) * beta(3) + strain_til(3) * beta(6) * beta(6) +
350 strain_til(4) * beta(0) * beta(6) + strain_til(5) * beta(3) * beta(6);
351 strain(1) = strain_til(0) * beta(1) * beta(1) + strain_til(1) * beta(4) * beta(4) +
352 strain_til(2) * beta(1) * beta(4) + strain_til(3) * beta(7) * beta(7) +
353 strain_til(4) * beta(1) * beta(7) + strain_til(5) * beta(4) * beta(7);
354 strain(2) = strain_til(0) * 2.0 * beta(0) * beta(1) + strain_til(1) * 2.0 * beta(3) * beta(4) +
355 strain_til(2) * (beta(1) * beta(3) + beta(0) * beta(4)) + strain_til(3) * 2.0 * beta(6) * beta(7) +
356 strain_til(4) * (beta(1) * beta(6) + beta(0) * beta(7)) +
357 strain_til(5) * (beta(4) * beta(6) + beta(3) * beta(7));
358 strain(3) = strain_til(0) * beta(2) * beta(2) + strain_til(1) * beta(5) * beta(5) +
359 strain_til(2) * beta(2) * beta(5) + strain_til(3) * beta(8) * beta(8) +
360 strain_til(4) * beta(2) * beta(8) + strain_til(5) * beta(5) * beta(8);
361 strain(4) = strain_til(0) * 2.0 * beta(0) * beta(2) + strain_til(1) * 2.0 * beta(3) * beta(5) +
362 strain_til(2) * (beta(2) * beta(3) + beta(0) * beta(5)) + strain_til(3) * 2.0 * beta(6) * beta(8) +
363 strain_til(4) * (beta(2) * beta(6) + beta(0) * beta(8)) +
364 strain_til(5) * (beta(5) * beta(6) + beta(3) * beta(8));
365 strain(5) = strain_til(0) * 2.0 * beta(1) * beta(2) + strain_til(1) * 2.0 * beta(4) * beta(5) +
366 strain_til(2) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3) * 2.0 * beta(7) * beta(8) +
367 strain_til(4) * (beta(2) * beta(7) + beta(1) * beta(8)) +
368 strain_til(5) * (beta(5) * beta(7) + beta(4) * beta(8));
369
370 // Straint derivative component
371
372 ChMatrixNM<double, 6, 24> strainD_til;
373 strainD_til.row(0) = Nx * (*d) * Sx;
374 strainD_til.row(1) = Ny * (*d) * Sy;
375 strainD_til.row(2) = Nx * (*d) * Sy + Ny * (*d) * Sx;
376 //== Compatible strain (No ANS)==//
377 strainD_til.row(3) = Nz * (*d) * Sz;
378 strainD_til.row(4) = Nx * (*d) * Sz + Nz * (*d) * Sx;
379 strainD_til.row(5) = Ny * (*d) * Sz + Nz * (*d) * Sy;
380
381 // For orthotropic material
382 for (int ii = 0; ii < 24; ii++) {
383 strainD(0, ii) = strainD_til(0, ii) * beta(0) * beta(0) + strainD_til(1, ii) * beta(3) * beta(3) +
384 strainD_til(2, ii) * beta(0) * beta(3) + strainD_til(3, ii) * beta(6) * beta(6) +
385 strainD_til(4, ii) * beta(0) * beta(6) + strainD_til(5, ii) * beta(3) * beta(6);
386 strainD(1, ii) = strainD_til(0, ii) * beta(1) * beta(1) + strainD_til(1, ii) * beta(4) * beta(4) +
387 strainD_til(2, ii) * beta(1) * beta(4) + strainD_til(3, ii) * beta(7) * beta(7) +
388 strainD_til(4, ii) * beta(1) * beta(7) + strainD_til(5, ii) * beta(4) * beta(7);
389 strainD(2, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(1) + strainD_til(1, ii) * 2.0 * beta(3) * beta(4) +
390 strainD_til(2, ii) * (beta(1) * beta(3) + beta(0) * beta(4)) +
391 strainD_til(3, ii) * 2.0 * beta(6) * beta(7) +
392 strainD_til(4, ii) * (beta(1) * beta(6) + beta(0) * beta(7)) +
393 strainD_til(5, ii) * (beta(4) * beta(6) + beta(3) * beta(7));
394 strainD(3, ii) = strainD_til(0, ii) * beta(2) * beta(2) + strainD_til(1, ii) * beta(5) * beta(5) +
395 strainD_til(2, ii) * beta(2) * beta(5) + strainD_til(3, ii) * beta(8) * beta(8) +
396 strainD_til(4, ii) * beta(2) * beta(8) + strainD_til(5) * beta(5) * beta(8);
397 strainD(4, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(2) + strainD_til(1, ii) * 2.0 * beta(3) * beta(5) +
398 strainD_til(2, ii) * (beta(2) * beta(3) + beta(0) * beta(5)) +
399 strainD_til(3, ii) * 2.0 * beta(6) * beta(8) +
400 strainD_til(4, ii) * (beta(2) * beta(6) + beta(0) * beta(8)) +
401 strainD_til(5, ii) * (beta(5) * beta(6) + beta(3) * beta(8));
402 strainD(5, ii) = strainD_til(0, ii) * 2.0 * beta(1) * beta(2) + strainD_til(1, ii) * 2.0 * beta(4) * beta(5) +
403 strainD_til(2, ii) * (beta(2) * beta(4) + beta(1) * beta(5)) +
404 strainD_til(3, ii) * 2.0 * beta(7) * beta(8) +
405 strainD_til(4, ii) * (beta(2) * beta(7) + beta(1) * beta(8)) +
406 strainD_til(5, ii) * (beta(5) * beta(7) + beta(4) * beta(8));
407 }
408
409 // Gd (8x24) calculation
410 for (int ii = 0; ii < 8; ii++) {
411 Gd(0, 3 * (ii) + 0) = j0(0, 0) * Nx(ii) + j0(1, 0) * Ny(ii) + j0(2, 0) * Nz(ii);
412 Gd(1, 3 * (ii) + 1) = j0(0, 0) * Nx(ii) + j0(1, 0) * Ny(ii) + j0(2, 0) * Nz(ii);
413 Gd(2, 3 * (ii) + 2) = j0(0, 0) * Nx(ii) + j0(1, 0) * Ny(ii) + j0(2, 0) * Nz(ii);
414
415 Gd(3, 3 * (ii) + 0) = j0(0, 1) * Nx(ii) + j0(1, 1) * Ny(ii) + j0(2, 1) * Nz(ii);
416 Gd(4, 3 * (ii) + 1) = j0(0, 1) * Nx(ii) + j0(1, 1) * Ny(ii) + j0(2, 1) * Nz(ii);
417 Gd(5, 3 * (ii) + 2) = j0(0, 1) * Nx(ii) + j0(1, 1) * Ny(ii) + j0(2, 1) * Nz(ii);
418
419 Gd(6, 3 * (ii) + 0) = j0(0, 2) * Nx(ii) + j0(1, 2) * Ny(ii) + j0(2, 2) * Nz(ii);
420 Gd(7, 3 * (ii) + 1) = j0(0, 2) * Nx(ii) + j0(1, 2) * Ny(ii) + j0(2, 2) * Nz(ii);
421 Gd(8, 3 * (ii) + 2) = j0(0, 2) * Nx(ii) + j0(1, 2) * Ny(ii) + j0(2, 2) * Nz(ii);
422 }
423
424 // Enhanced Assumed Strain 2nd
425 strain += strain_EAS;
426
427 ChMatrixNM<double, 9, 6> temp56;
428 ChVectorN<double, 9> HE1;
429 ChMatrixNM<double, 9, 24> GDEPSP;
430 ChMatrixNM<double, 9, 9> KALPHA;
431
432 // If Mooney-Rivlin Material is selected -> Calculates internal forces and their Jacobian accordingly (new E_eps)
433 if (element->m_isMooney) {
434 ChMatrix33<> CG; // CG: Right Cauchy-Green deformation tensor C=trans(F)*F
435 ChMatrix33<> INVCG; // INVCG: Inverse of Right Cauchy-Green deformation tensor C=trans(F)*F
436 ChMatrix33<> I1PC; // Stress tensor from first term of Mooney-Rivlin strain energy
437 ChMatrix33<> I2PC; // Stress tensor from second term of Mooney-Rivlin strain energy
438 ChMatrix33<> JPC; // Stress tensor from penalty term to ensure incompressibility
439 ChMatrix33<> STR; // Stress tensor from strain energy (with penalty for incompressibility CCOM3)
440
441 // Same quantities for the numerical calculation of the Jacobian of Mooney-Rivlin forces
442 ChMatrix33<> CGN; // CG: Right Cauchy-Green deformation tensor C=trans(F)*F
443 ChMatrix33<> INVCGN; // INVCG: Inverse of Right Cauchy-Green deformation tensor C=trans(F)*F
444 ChMatrix33<> I1PCN; // Stress tensor from first term of Mooney-Rivlin strain energy
445 ChMatrix33<> I2PCN; // Stress tensor from second term of Mooney-Rivlin strain energy
446 ChMatrix33<> JPCN; // Stress tensor from penalty term to ensure incompressibility
447 ChMatrix33<> STRN; // Stress tensor from strain energy (with penalty for incompressibility CCOM3)
448
449 ChVectorN<double, 6> strain_1;
450
451 // Right Cauchy - Green deformation tensor
452 CG(0, 0) = 2.0 * strain(0) + 1.0;
453 CG(1, 1) = 2.0 * strain(1) + 1.0;
454 CG(2, 2) = 2.0 * strain(3) + 1.0;
455 CG(1, 0) = strain(2);
456 CG(0, 1) = CG(1, 0);
457 CG(2, 0) = strain(4);
458 CG(0, 2) = CG(2, 0);
459 CG(2, 1) = strain(5);
460 CG(1, 2) = CG(2, 1);
461
462 INVCG = CG.inverse();
463
464 // Calculation of invariants I1, I2, and I3 and its deviatoric counterparts I1bar, I2bar, and I3bar
465 double Deld = 0.000001;
466 // First invariant of Right Cauchy-Green deformation tensor
467 double I1 = CG(0, 0) + CG(1, 1) + CG(2, 2);
468 // Second invariant of Right Cauchy-Green deformation tensor
469 double I2 = 0.5 * (pow(I1, 2) - (pow(CG(0, 0), 2) + pow(CG(1, 0), 2) + pow(CG(2, 0), 2) + pow(CG(0, 1), 2) +
470 pow(CG(1, 1), 2) + pow(CG(2, 1), 2) + pow(CG(0, 2), 2) + pow(CG(1, 2), 2) +
471 pow(CG(2, 2), 2)));
472 // Third invariant of Right Cauchy-Green deformation tensor (must be very close to 1 for incompressible
473 // material)
474 double I3 = CG(0, 0) * CG(1, 1) * CG(2, 2) - CG(0, 0) * CG(1, 2) * CG(2, 1) + CG(0, 1) * CG(1, 2) * CG(2, 0) -
475 CG(0, 1) * CG(1, 0) * CG(2, 2) + CG(0, 2) * CG(1, 0) * CG(2, 1) - CG(2, 0) * CG(1, 1) * CG(0, 2);
476 double J = sqrt(I3);
477 // double CCOM1 = 551584.0; // C10 not 0.551584
478 // double CCOM2 = 137896.0; // C01 not 0.137896
479 double CCOM3 = 2.0 * (element->CCOM1 + element->CCOM2) / (1.0 - 2.0 * 0.49); // K:bulk modulus
480 double StockEPS;
481
482 /// Calculation of stress tensor STR term to term: I1PC, I2PC, and JPC.
483
484 // Stress tensor from first term of Mooney-Rivlin strain energy
485 I1PC = (ChMatrix33<>::Identity() - INVCG * (1.0 / 3.0 * I1)) * pow(I3, -1.0 / 3.0);
486 // Stress tensor from second term of Mooney-Rivlin strain energy
487 I2PC = (((ChMatrix33<>::Identity() * I1) - CG) - (INVCG * (2.0 / 3.0) * I2)) * pow(I3, -2.0 / 3.0);
488 // Stress tensor from penalty for incompressibility
489 JPC = INVCG * (J / 2.0);
490 // Definition of stress tensor from strain energy (including penalty for incompressibility CCOM3)
491 STR = I1PC * (element->CCOM1 * 2.0) + I2PC * (element->CCOM2 * 2.0) + JPC * (CCOM3 * (J - 1.0) * 2.0);
492
493 // Put the stress in vector form
494 ChVectorN<double, 6> TEMP5;
495 TEMP5(0) = STR(0, 0);
496 TEMP5(1) = STR(1, 1);
497 TEMP5(2) = STR(0, 1);
498 TEMP5(3) = STR(2, 2);
499 TEMP5(4) = STR(0, 2);
500 TEMP5(5) = STR(1, 2);
501
502 E_eps.setZero();
503
504 // Compatible plus enhanced assumed strain
505 strain_1 = strain;
506
507 // Loop to obtain our Mooney-Rivling E_eps (tangential matrix of elastic coefficients)
508 // E_eps is necessary for obtaining the Jacobian of MR internal forces
509 ChVectorN<double, 6> TEMP5N;
510 for (int JJJ = 0; JJJ < 6; JJJ++) {
511 StockEPS = strain_1(JJJ);
512 strain_1(JJJ) = StockEPS + Deld;
513 CGN(0, 0) = 2.0 * strain_1(0) + 1.0;
514 CGN(1, 1) = 2.0 * strain_1(1) + 1.0;
515 CGN(2, 2) = 2.0 * strain_1(3) + 1.0;
516 CGN(1, 0) = strain_1(2);
517 CGN(0, 1) = CGN(1, 0);
518 CGN(2, 0) = strain_1(4);
519 CGN(0, 2) = CGN(2, 0);
520 CGN(2, 1) = strain_1(5);
521 CGN(1, 2) = CGN(2, 1);
522 INVCGN = CGN.inverse();
523
524 I1 = CGN(0, 0) + CGN(1, 1) + CGN(2, 2);
525 I2 = 0.5 * (pow(I1, 2) - (pow(CGN(0, 0), 2) + pow(CGN(1, 0), 2) + pow(CGN(2, 0), 2) + pow(CGN(0, 1), 2) +
526 pow(CGN(1, 1), 2) + pow(CGN(2, 1), 2) + pow(CGN(0, 2), 2) + pow(CGN(1, 2), 2) +
527 pow(CGN(2, 2), 2)));
528 I3 = CGN(0, 0) * CGN(1, 1) * CGN(2, 2) - CGN(0, 0) * CGN(1, 2) * CGN(2, 1) +
529 CGN(0, 1) * CGN(1, 2) * CGN(2, 0) - CGN(0, 1) * CGN(1, 0) * CGN(2, 2) +
530 CGN(0, 2) * CGN(1, 0) * CGN(2, 1) - CGN(2, 0) * CGN(1, 1) * CGN(0, 2);
531 J = sqrt(I3);
532 I1PCN = (ChMatrix33<>::Identity() - INVCGN * (1.0 / 3.0 * I1)) * pow(I3, -1.0 / 3.0);
533 I2PCN = (((ChMatrix33<>::Identity() * I1) - CGN) - (INVCGN * (2.0 / 3.0) * I2)) * pow(I3, -2.0 / 3.0);
534 JPCN = INVCGN * (J / 2.0);
535 STRN = I1PCN * (element->CCOM1 * 2.0) + I2PCN * (element->CCOM2 * 2.0) + JPCN * (CCOM3 * (J - 1.0) * 2.0);
536 TEMP5N(0) = STRN(0, 0);
537 TEMP5N(1) = STRN(1, 1);
538 TEMP5N(2) = STRN(0, 1);
539 TEMP5N(3) = STRN(2, 2);
540 TEMP5N(4) = STRN(0, 2);
541 TEMP5N(5) = STRN(1, 2);
542 strain_1(JJJ) = StockEPS;
543 E_eps(JJJ, 0) = (TEMP5N(0) - TEMP5(0)) / Deld;
544 E_eps(JJJ, 1) = (TEMP5N(1) - TEMP5(1)) / Deld;
545 E_eps(JJJ, 2) = (TEMP5N(2) - TEMP5(2)) / Deld;
546 E_eps(JJJ, 3) = (TEMP5N(3) - TEMP5(3)) / Deld;
547 E_eps(JJJ, 4) = (TEMP5N(4) - TEMP5(4)) / Deld;
548 E_eps(JJJ, 5) = (TEMP5N(5) - TEMP5(5)) / Deld;
549 }
550 // Add internal forces to Fint and HE1 for Mooney-Rivlin
551 temp56 = G.transpose() * E_eps;
552 Fint = strainD.transpose() * TEMP5;
553 Fint *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
554 HE1 = G.transpose() * TEMP5;
555 HE1 *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
556 Sigm(0, 0) = TEMP5(0);
557 Sigm(1, 1) = TEMP5(0);
558 Sigm(2, 2) = TEMP5(0);
559 Sigm(0, 3) = TEMP5(2);
560 Sigm(1, 4) = TEMP5(2);
561 Sigm(2, 5) = TEMP5(2);
562 Sigm(0, 6) = TEMP5(4);
563 Sigm(1, 7) = TEMP5(4);
564 Sigm(2, 8) = TEMP5(4);
565 Sigm(3, 0) = TEMP5(2);
566 Sigm(4, 1) = TEMP5(2);
567 Sigm(5, 2) = TEMP5(2);
568 Sigm(3, 3) = TEMP5(1);
569 Sigm(4, 4) = TEMP5(1);
570 Sigm(5, 5) = TEMP5(1);
571 Sigm(3, 6) = TEMP5(5);
572 Sigm(4, 7) = TEMP5(5);
573 Sigm(5, 8) = TEMP5(5);
574 Sigm(6, 0) = TEMP5(4);
575 Sigm(7, 1) = TEMP5(4);
576 Sigm(8, 2) = TEMP5(4);
577 Sigm(6, 3) = TEMP5(5);
578 Sigm(7, 4) = TEMP5(5);
579 Sigm(8, 5) = TEMP5(5);
580 Sigm(6, 6) = TEMP5(3);
581 Sigm(7, 7) = TEMP5(3);
582 Sigm(8, 8) = TEMP5(3);
583 } else {
584 /// Stress tensor calculation
585 stress = E_eps * strain;
586 Sigm(0, 0) = stress(0);
587 Sigm(0, 3) = stress(2);
588 Sigm(0, 6) = stress(4);
589 Sigm(1, 1) = stress(0);
590 Sigm(1, 4) = stress(2);
591 Sigm(1, 7) = stress(4);
592 Sigm(2, 2) = stress(0);
593 Sigm(2, 5) = stress(2);
594 Sigm(2, 8) = stress(4);
595 // XX //XY //XZ
596 Sigm(3, 0) = stress(2);
597 Sigm(3, 3) = stress(1);
598 Sigm(3, 6) = stress(5);
599 Sigm(4, 1) = stress(2);
600 Sigm(4, 4) = stress(1);
601 Sigm(4, 7) = stress(5);
602 Sigm(5, 2) = stress(2);
603 Sigm(5, 5) = stress(1);
604 Sigm(5, 8) = stress(5);
605 // XY //YY //YZ
606 Sigm(6, 0) = stress(4);
607 Sigm(6, 3) = stress(5);
608 Sigm(6, 6) = stress(3);
609 Sigm(7, 1) = stress(4);
610 Sigm(7, 4) = stress(5);
611 Sigm(7, 7) = stress(3);
612 Sigm(8, 2) = stress(4);
613 Sigm(8, 5) = stress(5);
614 Sigm(8, 8) = stress(3);
615 // XZ //YZ //ZZ
616 // Add internal forces to Fint and HE1 for linear elastic material
617 // temp56 is actually 9x6 for the brick element (9 EAS internal parameters)
618 temp56 = G.transpose() * E_eps;
619 // Add generalized internal force
620 Fint = strainD.transpose() * E_eps * strain;
621 Fint *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
622 // Add EAS internal force (vector of 9 components for each element)
623 HE1 = temp56 * strain;
624 HE1 *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
625 } // end of if(isMooney==1)
626
627 // Internal force (linear isotropic or Mooney-Rivlin) Jacobian calculation
628 // First term for Jacobian matrix
629 // Second term for Jacobian matrix
630 JAC11 = strainD.transpose() * E_eps * strainD + Gd.transpose() * Sigm * Gd;
631 // Final expression for the Jacobian
632 JAC11 *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
633 // Jacobian of EAS forces w.r.t. element coordinates
634 double factor_g =
635 detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
636 GDEPSP = factor_g * temp56 * strainD;
637 // Jacobian of EAS forces (w.r.t. EAS internal parameters)
638 double factor_k =
639 detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
640 KALPHA = factor_k * temp56 * G;
641
642 ChVectorN<double, 216> GDEPSPVec;
643 ChVectorN<double, 81> KALPHAVec;
644 ChVectorN<double, 576> JACVec;
645
646 for (int i = 0; i < 9; i++) {
647 for (int j = 0; j < 24; j++) {
648 GDEPSPVec(i * 24 + j) = GDEPSP(i, j);
649 }
650 }
651
652 // GDEPSP = GDEPSPvec;
653 for (int i = 0; i < 9; i++) {
654 for (int j = 0; j < 9; j++) {
655 KALPHAVec(i * 9 + j) = KALPHA(i, j);
656 }
657 }
658
659 // KALPHAVec = KALPHA;
660 for (int i = 0; i < 24; i++) {
661 for (int j = 0; j < 24; j++) {
662 JACVec(i * 24 + j) = JAC11(i, j);
663 }
664 }
665
666 // JACVec = JAC11;
667 result.segment(0, 24) = Fint;
668 result.segment(24, 9) = HE1;
669 result.segment(33, 216) = GDEPSPVec;
670 result.segment(249, 81) = KALPHAVec;
671 result.segment(330, 576) = JACVec;
672 }
673
674 // -----------------------------------------------------------------------------
675
676 class Brick_ForceNumerical : public ChIntegrable3D<ChVectorN<double, 330>> {
677 public:
678 Brick_ForceNumerical(ChMatrixNM<double, 8, 3>* d_,
679 ChMatrixNM<double, 8, 3>* d0_,
680 ChElementHexaANCF_3813* element_,
681 ChMatrixNM<double, 6, 6>* T0_,
682 double* detJ0C_,
683 ChVectorN<double, 9>* alpha_eas_);
684
685 Brick_ForceNumerical(ChMatrixNM<double, 8, 3>* d_,
686 ChMatrixNM<double, 8, 3>* d0_,
687 ChElementHexaANCF_3813* element_,
688 ChMatrixNM<double, 6, 6>* T0_,
689 double* detJ0C_,
690 ChVectorN<double, 9>* alpha_eas_,
691 double* E_,
692 double* v_);
~Brick_ForceNumerical()693 ~Brick_ForceNumerical() {}
694
695 private:
696 ChElementHexaANCF_3813* element;
697 // Pointers used for external values
698 ChMatrixNM<double, 8, 3>* d; // Pointer to a matrix containing the element coordinates
699 ChMatrixNM<double, 8, 3>* d0; // Pointer to a matrix containing the element initial coordinates
700 ChMatrixNM<double, 6, 6>* T0; // Pointer to transformation matrix for Enhanced Assumed Strain (EAS)
701 ChVectorN<double, 9>* alpha_eas; // Pointer to the 9 internal parameters for EAS
702 double* detJ0C; // Pointer to determinant of the initial Jacobian at the element center
703 double* E; // Pointer to Young modulus
704 double* v; // Pointer to Poisson ratio
705 ChVectorN<double, 24> Fint; // Generalized internal (elastic) force vector
706 ChMatrixNM<double, 6, 6> E_eps; // Matrix of elastic coefficients (features orthotropy)
707 ChMatrixNM<double, 3, 24> Sx; // Sparse shape function matrix, X derivative
708 ChMatrixNM<double, 3, 24> Sy; // Sparse shape function matrix, Y derivative
709 ChMatrixNM<double, 3, 24> Sz; // Sparse shape function matrix, Z derivative
710 ChElementHexaANCF_3813::ShapeVector Nx; // Dense shape function vector, X derivative
711 ChElementHexaANCF_3813::ShapeVector Ny; // Dense shape function vector, Y derivative
712 ChElementHexaANCF_3813::ShapeVector Nz; // Dense shape function vector, Z derivative
713 ChMatrixNM<double, 6, 24> strainD; // Derivative of the strains w.r.t. the coordinates. Includes orthotropy
714 ChVectorN<double, 6> strain; // Vector of strains
715 double detJ0; // Determinant of the initial position vector gradient matrix
716 // EAS
717 ChMatrixNM<double, 6, 9> M; // Shape function matrix for Enhanced Assumed Strain
718 ChMatrixNM<double, 6, 9> G; // Matrix G interpolates the internal parameters of EAS
719 ChVectorN<double, 6> strain_EAS; // Enhanced assumed strain vector
720
721 // Gaussian integration to calculate internal forces and EAS matrices
722 virtual void Evaluate(ChVectorN<double, 330>& result, const double x, const double y, const double z) override;
723 };
724
Brick_ForceNumerical(ChMatrixNM<double,8,3> * d_,ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_,ChMatrixNM<double,6,6> * T0_,double * detJ0C_,ChVectorN<double,9> * alpha_eas_)725 Brick_ForceNumerical::Brick_ForceNumerical(ChMatrixNM<double, 8, 3>* d_,
726 ChMatrixNM<double, 8, 3>* d0_,
727 ChElementHexaANCF_3813* element_,
728 ChMatrixNM<double, 6, 6>* T0_,
729 double* detJ0C_,
730 ChVectorN<double, 9>* alpha_eas_)
731 : element(element_), d(d_), d0(d0_), T0(T0_), alpha_eas(alpha_eas_), detJ0C(detJ0C_) {
732 E_eps.setZero();
733
734 Sx.setZero();
735 Sy.setZero();
736 Sz.setZero();
737 }
738
Brick_ForceNumerical(ChMatrixNM<double,8,3> * d_,ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_,ChMatrixNM<double,6,6> * T0_,double * detJ0C_,ChVectorN<double,9> * alpha_eas_,double * E_,double * v_)739 Brick_ForceNumerical::Brick_ForceNumerical(ChMatrixNM<double, 8, 3>* d_,
740 ChMatrixNM<double, 8, 3>* d0_,
741 ChElementHexaANCF_3813* element_,
742 ChMatrixNM<double, 6, 6>* T0_,
743 double* detJ0C_,
744 ChVectorN<double, 9>* alpha_eas_,
745 double* E_,
746 double* v_)
747 : element(element_), d(d_), d0(d0_), T0(T0_), alpha_eas(alpha_eas_), detJ0C(detJ0C_), E(E_), v(v_) {
748 E_eps.setZero();
749
750 Sx.setZero();
751 Sy.setZero();
752 Sz.setZero();
753 }
754
Evaluate(ChVectorN<double,330> & result,const double x,const double y,const double z)755 void Brick_ForceNumerical::Evaluate(ChVectorN<double, 330>& result, const double x, const double y, const double z) {
756 element->ShapeFunctionsDerivativeX(Nx, x, y, z);
757 element->ShapeFunctionsDerivativeY(Ny, x, y, z);
758 element->ShapeFunctionsDerivativeZ(Nz, x, y, z);
759 element->Basis_M(M, x, y, z); // EAS
760
761 if (!element->m_isMooney) { // m_isMooney == false means linear elastic material
762 double DD = (*E) * (1.0 - (*v)) / ((1.0 + (*v)) * (1.0 - 2.0 * (*v)));
763 E_eps.fillDiagonal(1.0);
764 E_eps(0, 1) = (*v) / (1.0 - (*v));
765 E_eps(0, 3) = (*v) / (1.0 - (*v));
766 E_eps(1, 0) = (*v) / (1.0 - (*v));
767 E_eps(1, 3) = (*v) / (1.0 - (*v));
768 E_eps(2, 2) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
769 E_eps(3, 0) = (*v) / (1.0 - (*v));
770 E_eps(3, 1) = (*v) / (1.0 - (*v));
771 E_eps(4, 4) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
772 E_eps(5, 5) = (1.0 - 2.0 * (*v)) / (2.0 * (1.0 - (*v)));
773 E_eps *= DD;
774 }
775 // Expand shape functions Sx, Sy, Sz
776 // Sx=[Nx1*eye(3) Nx2*eye(3) Nx3*eye(3) Nx4*eye(3) Nx5*eye(3) Nx6*eye(3) Nx7*eye(3) Nx8*eye(3)]
777
778 for (int i = 0; i < 8; i++) {
779 Sx(0, 3 * i + 0) = Nx(i);
780 Sx(1, 3 * i + 1) = Nx(i);
781 Sx(2, 3 * i + 2) = Nx(i);
782 }
783
784 for (int i = 0; i < 8; i++) {
785 Sy(0, 3 * i + 0) = Ny(i);
786 Sy(1, 3 * i + 1) = Ny(i);
787 Sy(2, 3 * i + 2) = Ny(i);
788 }
789
790 for (int i = 0; i < 8; i++) {
791 Sz(0, 3 * i + 0) = Nz(i);
792 Sz(1, 3 * i + 1) = Nz(i);
793 Sz(2, 3 * i + 2) = Nz(i);
794 }
795
796 //==EAS and Initial Shape==//
797 ChMatrixNM<double, 3, 3> rd0;
798 rd0.col(0) = (*d0).transpose() * Nx.transpose();
799 rd0.col(1) = (*d0).transpose() * Ny.transpose();
800 rd0.col(2) = (*d0).transpose() * Nz.transpose();
801 detJ0 = rd0.determinant();
802
803 //////////////////////////////////////////////////////////////
804 //// Transformation : Orthogonal transformation (A and J) ////
805 //////////////////////////////////////////////////////////////
806 ChVector<double> G1;
807 ChVector<double> G2;
808 ChVector<double> G3;
809 ChVector<double> G1xG2;
810 G1[0] = rd0(0, 0);
811 G2[0] = rd0(0, 1);
812 G3[0] = rd0(0, 2);
813 G1[1] = rd0(1, 0);
814 G2[1] = rd0(1, 1);
815 G3[1] = rd0(1, 2);
816 G1[2] = rd0(2, 0);
817 G2[2] = rd0(2, 1);
818 G3[2] = rd0(2, 2);
819 G1xG2.Cross(G1, G2);
820
821 ////Tangent Frame
822 ChVector<> A1 = G1 / sqrt(G1[0] * G1[0] + G1[1] * G1[1] + G1[2] * G1[2]);
823 ChVector<> A3 = G1xG2 / sqrt(G1xG2[0] * G1xG2[0] + G1xG2[1] * G1xG2[1] + G1xG2[2] * G1xG2[2]);
824 ChVector<> A2 = A3.Cross(A1);
825
826 ////Direction for orthotropic material//
827 double theta = 0.0;
828 ChVector<> AA1 = A1 * cos(theta) + A2 * sin(theta);
829 ChVector<> AA2 = -A1 * sin(theta) + A2 * cos(theta);
830 ChVector<> AA3 = A3;
831
832 ////Beta
833 ChMatrixNM<double, 3, 3> j0 = rd0.inverse();
834 ChVector<double> j01;
835 ChVector<double> j02;
836 ChVector<double> j03;
837 j01[0] = j0(0, 0);
838 j02[0] = j0(1, 0);
839 j03[0] = j0(2, 0);
840 j01[1] = j0(0, 1);
841 j02[1] = j0(1, 1);
842 j03[1] = j0(2, 1);
843 j01[2] = j0(0, 2);
844 j02[2] = j0(1, 2);
845 j03[2] = j0(2, 2);
846 ChVectorN<double, 9> beta;
847 beta(0) = Vdot(AA1, j01);
848 beta(1) = Vdot(AA2, j01);
849 beta(2) = Vdot(AA3, j01);
850 beta(3) = Vdot(AA1, j02);
851 beta(4) = Vdot(AA2, j02);
852 beta(5) = Vdot(AA3, j02);
853 beta(6) = Vdot(AA1, j03);
854 beta(7) = Vdot(AA2, j03);
855 beta(8) = Vdot(AA3, j03);
856
857 //////////////////////////////////////////////////
858 //// Enhanced Assumed Strain /////////////////////
859 //////////////////////////////////////////////////
860 G = (*T0) * M * ((*detJ0C) / detJ0);
861 strain_EAS = G * (*alpha_eas);
862
863 ChMatrixNM<double, 8, 8> d_d = (*d) * (*d).transpose();
864 ChVectorN<double, 8> ddNx = d_d * Nx.transpose();
865 ChVectorN<double, 8> ddNy = d_d * Ny.transpose();
866 ChVectorN<double, 8> ddNz = d_d * Nz.transpose();
867
868 ChMatrixNM<double, 8, 8> d0_d0 = (*d0) * (*d0).transpose();
869 ChVectorN<double, 8> d0d0Nx = d0_d0 * Nx.transpose();
870 ChVectorN<double, 8> d0d0Ny = d0_d0 * Ny.transpose();
871 ChVectorN<double, 8> d0d0Nz = d0_d0 * Nz.transpose();
872
873 ///////////////////////////
874 /// Strain component //////
875 ///////////////////////////
876 ChVectorN<double, 6> strain_til;
877 strain_til(0) = 0.5 * Nx.dot(ddNx - d0d0Nx);
878 strain_til(1) = 0.5 * Ny.dot(ddNy - d0d0Ny);
879 strain_til(2) = Nx.dot(ddNy - d0d0Ny);
880 //== Compatible strain (No ANS) ==//
881 strain_til(3) = 0.5 * Nz.dot(ddNz - d0d0Nz);
882 strain_til(4) = Nx.dot(ddNz - d0d0Nz);
883 strain_til(5) = Ny.dot(ddNz - d0d0Nz);
884
885 //// For orthotropic material ///
886 strain(0) = strain_til(0) * beta(0) * beta(0) + strain_til(1) * beta(3) * beta(3) +
887 strain_til(2) * beta(0) * beta(3) + strain_til(3) * beta(6) * beta(6) +
888 strain_til(4) * beta(0) * beta(6) + strain_til(5) * beta(3) * beta(6);
889 strain(1) = strain_til(0) * beta(1) * beta(1) + strain_til(1) * beta(4) * beta(4) +
890 strain_til(2) * beta(1) * beta(4) + strain_til(3) * beta(7) * beta(7) +
891 strain_til(4) * beta(1) * beta(7) + strain_til(5) * beta(4) * beta(7);
892 strain(2) = strain_til(0) * 2.0 * beta(0) * beta(1) + strain_til(1) * 2.0 * beta(3) * beta(4) +
893 strain_til(2) * (beta(1) * beta(3) + beta(0) * beta(4)) + strain_til(3) * 2.0 * beta(6) * beta(7) +
894 strain_til(4) * (beta(1) * beta(6) + beta(0) * beta(7)) +
895 strain_til(5) * (beta(4) * beta(6) + beta(3) * beta(7));
896 strain(3) = strain_til(0) * beta(2) * beta(2) + strain_til(1) * beta(5) * beta(5) +
897 strain_til(2) * beta(2) * beta(5) + strain_til(3) * beta(8) * beta(8) +
898 strain_til(4) * beta(2) * beta(8) + strain_til(5) * beta(5) * beta(8);
899 strain(4) = strain_til(0) * 2.0 * beta(0) * beta(2) + strain_til(1) * 2.0 * beta(3) * beta(5) +
900 strain_til(2) * (beta(2) * beta(3) + beta(0) * beta(5)) + strain_til(3) * 2.0 * beta(6) * beta(8) +
901 strain_til(4) * (beta(2) * beta(6) + beta(0) * beta(8)) +
902 strain_til(5) * (beta(5) * beta(6) + beta(3) * beta(8));
903 strain(5) = strain_til(0) * 2.0 * beta(1) * beta(2) + strain_til(1) * 2.0 * beta(4) * beta(5) +
904 strain_til(2) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3) * 2.0 * beta(7) * beta(8) +
905 strain_til(4) * (beta(2) * beta(7) + beta(1) * beta(8)) +
906 strain_til(5) * (beta(5) * beta(7) + beta(4) * beta(8));
907
908 ////////////////////////////////////
909 /// Straint derivative component ///
910 ////////////////////////////////////
911 ChMatrixNM<double, 6, 24> strainD_til;
912 strainD_til.row(0) = Nx * (*d) * Sx;
913 strainD_til.row(1) = Ny * (*d) * Sy;
914 strainD_til.row(2) = Nx * (*d) * Sy + Ny * (*d) * Sx;
915 //== Compatible strain (No ANS)==//
916 strainD_til.row(3) = Nz * (*d) * Sz;
917 strainD_til.row(4) = Nx * (*d) * Sz + Nz * (*d) * Sx;
918 strainD_til.row(5) = Ny * (*d) * Sz + Nz * (*d) * Sy;
919
920 //// For orthotropic material ///
921 for (int ii = 0; ii < 24; ii++) {
922 strainD(0, ii) = strainD_til(0, ii) * beta(0) * beta(0) + strainD_til(1, ii) * beta(3) * beta(3) +
923 strainD_til(2, ii) * beta(0) * beta(3) + strainD_til(3, ii) * beta(6) * beta(6) +
924 strainD_til(4, ii) * beta(0) * beta(6) + strainD_til(5, ii) * beta(3) * beta(6);
925 strainD(1, ii) = strainD_til(0, ii) * beta(1) * beta(1) + strainD_til(1, ii) * beta(4) * beta(4) +
926 strainD_til(2, ii) * beta(1) * beta(4) + strainD_til(3, ii) * beta(7) * beta(7) +
927 strainD_til(4, ii) * beta(1) * beta(7) + strainD_til(5, ii) * beta(4) * beta(7);
928 strainD(2, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(1) + strainD_til(1, ii) * 2.0 * beta(3) * beta(4) +
929 strainD_til(2, ii) * (beta(1) * beta(3) + beta(0) * beta(4)) +
930 strainD_til(3, ii) * 2.0 * beta(6) * beta(7) +
931 strainD_til(4, ii) * (beta(1) * beta(6) + beta(0) * beta(7)) +
932 strainD_til(5, ii) * (beta(4) * beta(6) + beta(3) * beta(7));
933 strainD(3, ii) = strainD_til(0, ii) * beta(2) * beta(2) + strainD_til(1, ii) * beta(5) * beta(5) +
934 strainD_til(2, ii) * beta(2) * beta(5) + strainD_til(3, ii) * beta(8) * beta(8) +
935 strainD_til(4, ii) * beta(2) * beta(8) + strainD_til(5) * beta(5) * beta(8);
936 strainD(4, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(2) + strainD_til(1, ii) * 2.0 * beta(3) * beta(5) +
937 strainD_til(2, ii) * (beta(2) * beta(3) + beta(0) * beta(5)) +
938 strainD_til(3, ii) * 2.0 * beta(6) * beta(8) +
939 strainD_til(4, ii) * (beta(2) * beta(6) + beta(0) * beta(8)) +
940 strainD_til(5, ii) * (beta(5) * beta(6) + beta(3) * beta(8));
941 strainD(5, ii) = strainD_til(0, ii) * 2.0 * beta(1) * beta(2) + strainD_til(1, ii) * 2.0 * beta(4) * beta(5) +
942 strainD_til(2, ii) * (beta(2) * beta(4) + beta(1) * beta(5)) +
943 strainD_til(3, ii) * 2.0 * beta(7) * beta(8) +
944 strainD_til(4, ii) * (beta(2) * beta(7) + beta(1) * beta(8)) +
945 strainD_til(5, ii) * (beta(5) * beta(7) + beta(4) * beta(8));
946 }
947
948 ///////////////////////////////////
949 /// Enhanced Assumed Strain 2nd ///
950 ///////////////////////////////////
951 strain += strain_EAS; // same as EPS in FORTRAN
952
953 ChMatrixNM<double, 9, 6> temp56; // same as TEMP1 in FORTRAN
954 ChVectorN<double, 9> HE1;
955 ChMatrixNM<double, 9, 24> GDEPSP;
956 ChMatrixNM<double, 9, 9> KALPHA;
957
958 // m_isMooney == 1 use Iso_Nonlinear_Mooney-Rivlin Material (2-parameters=> 3 inputs)
959 if (element->m_isMooney == 1) {
960 ChMatrix33<> CG; // CG: Right Cauchy-Green tensor C=trans(F)*F
961 ChMatrix33<> INVCG;
962 ChMatrix33<> I1PC;
963 ChMatrix33<> I2PC;
964 ChMatrix33<> JPC;
965 ChMatrix33<> STR;
966
967 ChMatrix33<> CGN;
968 ChMatrix33<> INVCGN;
969 ChMatrix33<> I1PCN;
970 ChMatrix33<> I2PCN;
971 ChMatrix33<> JPCN;
972 ChMatrix33<> STRN;
973
974 ChVectorN<double, 6> strain_1;
975
976 CG(0, 0) = 2.0 * strain(0) + 1.0;
977 CG(1, 1) = 2.0 * strain(1) + 1.0;
978 CG(2, 2) = 2.0 * strain(3) + 1.0;
979 CG(1, 0) = strain(2);
980 CG(0, 1) = CG(1, 0);
981 CG(2, 0) = strain(4);
982 CG(0, 2) = CG(2, 0);
983 CG(2, 1) = strain(5);
984 CG(1, 2) = CG(2, 1);
985
986 INVCG = CG.inverse();
987
988 double Deld = 0.000001;
989 double I1 = CG(0, 0) + CG(1, 1) + CG(2, 2);
990 double I2 = 0.5 * (pow(I1, 2) - (pow(CG(0, 0), 2) + pow(CG(1, 0), 2) + pow(CG(2, 0), 2) + pow(CG(0, 1), 2) +
991 pow(CG(1, 1), 2) + pow(CG(2, 1), 2) + pow(CG(0, 2), 2) + pow(CG(1, 2), 2) +
992 pow(CG(2, 2), 2)));
993 double I3 = CG(0, 0) * CG(1, 1) * CG(2, 2) - CG(0, 0) * CG(1, 2) * CG(2, 1) + CG(0, 1) * CG(1, 2) * CG(2, 0) -
994 CG(0, 1) * CG(1, 0) * CG(2, 2) + CG(0, 2) * CG(1, 0) * CG(2, 1) - CG(2, 0) * CG(1, 1) * CG(0, 2);
995 double J = sqrt(I3);
996 // double CCOM1 = 551584.0; // C10 not 0.551584
997 // double CCOM2 = 137896.0; // C01 not 0.137896
998 double CCOM3 = 2.0 * (element->CCOM1 + element->CCOM2) / (1.0 - 2.0 * 0.49); // K:bulk modulus
999 double StockEPS;
1000
1001 I1PC = (ChMatrix33<>::Identity() - INVCG * (1.0 / 3.0 * I1)) * pow(I3, -1.0 / 3.0);
1002 I2PC = (((ChMatrix33<>::Identity() * I1) - CG) - (INVCG * (2.0 / 3.0) * I2)) * pow(I3, -2.0 / 3.0);
1003 JPC = INVCG * (J / 2.0);
1004
1005 STR = I1PC * (element->CCOM1 * 2.0) + I2PC * (element->CCOM2 * 2.0) + JPC * (CCOM3 * (J - 1.0) * 2.0);
1006
1007 ChVectorN<double, 6> TEMP5;
1008 TEMP5(0) = STR(0, 0);
1009 TEMP5(1) = STR(1, 1);
1010 TEMP5(2) = STR(0, 1);
1011 TEMP5(3) = STR(2, 2);
1012 TEMP5(4) = STR(0, 2);
1013 TEMP5(5) = STR(1, 2);
1014
1015 E_eps.setZero();
1016
1017 ChVectorN<double, 6> TEMP5N;
1018 strain_1 = strain;
1019 for (int JJJ = 0; JJJ < 6; JJJ++) {
1020 StockEPS = strain_1(JJJ);
1021 strain_1(JJJ) = StockEPS + Deld;
1022 CGN(0, 0) = 2.0 * strain_1(0) + 1.0;
1023 CGN(1, 1) = 2.0 * strain_1(1) + 1.0;
1024 CGN(2, 2) = 2.0 * strain_1(3) + 1.0;
1025 CGN(1, 0) = strain_1(2);
1026 CGN(0, 1) = CGN(1, 0);
1027 CGN(2, 0) = strain_1(4);
1028 CGN(0, 2) = CGN(2, 0);
1029 CGN(2, 1) = strain_1(5);
1030 CGN(1, 2) = CGN(2, 1);
1031 INVCGN = CGN.inverse();
1032 I1 = CGN(0, 0) + CGN(1, 1) + CGN(2, 2);
1033 I2 = 0.5 * (pow(I1, 2) - (pow(CGN(0, 0), 2) + pow(CGN(1, 0), 2) + pow(CGN(2, 0), 2) + pow(CGN(0, 1), 2) +
1034 pow(CGN(1, 1), 2) + pow(CGN(2, 1), 2) + pow(CGN(0, 2), 2) + pow(CGN(1, 2), 2) +
1035 pow(CGN(2, 2), 2)));
1036 I3 = CGN(0, 0) * CGN(1, 1) * CGN(2, 2) - CGN(0, 0) * CGN(1, 2) * CGN(2, 1) +
1037 CGN(0, 1) * CGN(1, 2) * CGN(2, 0) - CGN(0, 1) * CGN(1, 0) * CGN(2, 2) +
1038 CGN(0, 2) * CGN(1, 0) * CGN(2, 1) - CGN(2, 0) * CGN(1, 1) * CGN(0, 2);
1039 J = sqrt(I3);
1040 I1PCN = (ChMatrix33<>::Identity() - INVCGN * (1.0 / 3.0 * I1)) * pow(I3, -1.0 / 3.0);
1041 I2PCN = (((ChMatrix33<>::Identity() * I1) - CGN) - (INVCGN * (2.0 / 3.0) * I2)) * pow(I3, -2.0 / 3.0);
1042 JPCN = INVCGN * (J / 2.0);
1043 STRN = I1PCN * (element->CCOM1 * 2.0) + I2PCN * (element->CCOM2 * 2.0) + JPCN * (CCOM3 * (J - 1.0) * 2.0);
1044 TEMP5N(0) = STRN(0, 0);
1045 TEMP5N(1) = STRN(1, 1);
1046 TEMP5N(2) = STRN(0, 1);
1047 TEMP5N(3) = STRN(2, 2);
1048 TEMP5N(4) = STRN(0, 2);
1049 TEMP5N(5) = STRN(1, 2);
1050 strain_1(JJJ) = StockEPS;
1051 E_eps(JJJ, 0) = (TEMP5N(0) - TEMP5(0)) / Deld;
1052 E_eps(JJJ, 1) = (TEMP5N(1) - TEMP5(1)) / Deld;
1053 E_eps(JJJ, 2) = (TEMP5N(2) - TEMP5(2)) / Deld;
1054 E_eps(JJJ, 3) = (TEMP5N(3) - TEMP5(3)) / Deld;
1055 E_eps(JJJ, 4) = (TEMP5N(4) - TEMP5(4)) / Deld;
1056 E_eps(JJJ, 5) = (TEMP5N(5) - TEMP5(5)) / Deld;
1057 }
1058 temp56 = G.transpose() * E_eps;
1059 Fint = strainD.transpose() * TEMP5;
1060 Fint *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1061 HE1 = G.transpose() * TEMP5;
1062 HE1 *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1063 } else {
1064 temp56 = G.transpose() * E_eps;
1065 Fint = strainD.transpose() * E_eps * strain;
1066 Fint *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1067 HE1 = temp56 * strain;
1068 HE1 *= detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1069 } // end of if(*flag_Mooney==1)
1070
1071 double factor_k =
1072 detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1073 KALPHA = factor_k * temp56 * G;
1074 double factor_g =
1075 detJ0 * (element->GetLengthX() / 2.0) * (element->GetLengthY() / 2.0) * (element->GetLengthZ() / 2.0);
1076 GDEPSP = factor_g * temp56 * strainD;
1077
1078 ChVectorN<double, 216> GDEPSPVec = Eigen::Map<ChVectorN<double, 216>>(GDEPSP.data(), 216);
1079 ChVectorN<double, 81> KALPHAVec = Eigen::Map<ChVectorN<double, 81>>(KALPHA.data(), 81);
1080
1081 result.segment(0, 24) = Fint;
1082 result.segment(24, 9) = HE1;
1083 result.segment(33, 216) = GDEPSPVec;
1084 result.segment(249, 81) = KALPHAVec;
1085 }
1086
1087 // -----------------------------------------------------------------------------
1088
ComputeInternalForces(ChVectorDynamic<> & Fi)1089 void ChElementHexaANCF_3813::ComputeInternalForces(ChVectorDynamic<>& Fi) {
1090 int ie = GetElemNum();
1091
1092 ChVector<> pA = m_nodes[0]->GetPos();
1093 ChVector<> pB = m_nodes[1]->GetPos();
1094 ChVector<> pC = m_nodes[2]->GetPos();
1095 ChVector<> pD = m_nodes[3]->GetPos();
1096 ChVector<> pE = m_nodes[4]->GetPos();
1097 ChVector<> pF = m_nodes[5]->GetPos();
1098 ChVector<> pG = m_nodes[6]->GetPos();
1099 ChVector<> pH = m_nodes[7]->GetPos();
1100
1101 ChMatrixNM<double, 8, 3> d;
1102 d(0, 0) = pA.x();
1103 d(0, 1) = pA.y();
1104 d(0, 2) = pA.z();
1105 d(1, 0) = pB.x();
1106 d(1, 1) = pB.y();
1107 d(1, 2) = pB.z();
1108 d(2, 0) = pC.x();
1109 d(2, 1) = pC.y();
1110 d(2, 2) = pC.z();
1111 d(3, 0) = pD.x();
1112 d(3, 1) = pD.y();
1113 d(3, 2) = pD.z();
1114 d(4, 0) = pE.x();
1115 d(4, 1) = pE.y();
1116 d(4, 2) = pE.z();
1117 d(5, 0) = pF.x();
1118 d(5, 1) = pF.y();
1119 d(5, 2) = pF.z();
1120 d(6, 0) = pG.x();
1121 d(6, 1) = pG.y();
1122 d(6, 2) = pG.z();
1123 d(7, 0) = pH.x();
1124 d(7, 1) = pH.y();
1125 d(7, 2) = pH.z();
1126
1127 double v = m_Material->Get_v();
1128 double E = m_Material->Get_E();
1129
1130 Fi.setZero();
1131
1132 /// If numerical differentiation is used, only the internal force and EAS stiffness
1133 /// will be calculated. If the numerical differentiation is not used, the jacobian
1134 /// will also be calculated.
1135 bool use_numerical_differentiation = false;
1136
1137 /// Internal force and EAS parameters are calculated for numerical differentiation.
1138 if (use_numerical_differentiation) {
1139 //////////////////////////////////////////////////////////////////////////////////////////////////////////
1140 ChVectorN<double, 330> TempIntegratedResult;
1141 ChVectorN<double, 24> Finternal;
1142
1143 ChMatrixNM<double, 6, 6> T0;
1144 ChVectorN<double, 9> HE;
1145 ChMatrixNM<double, 9, 24> GDEPSP;
1146 ChMatrixNM<double, 9, 9> KALPHA;
1147 ChMatrixNM<double, 9, 9> KALPHA1;
1148 ChVectorN<double, 9> ResidHE;
1149 double detJ0C;
1150 ChVectorN<double, 9> alpha_eas;
1151 ChVectorN<double, 9> renewed_alpha_eas;
1152 ChVectorN<double, 9> previous_alpha;
1153
1154 previous_alpha = m_stock_alpha_EAS;
1155 alpha_eas = previous_alpha;
1156 ResidHE.setZero();
1157 int count = 0;
1158 int fail = 1;
1159 /// Begin EAS loop
1160 while (fail == 1) {
1161 /// Update alpha EAS
1162 alpha_eas = alpha_eas - ResidHE;
1163 renewed_alpha_eas = alpha_eas;
1164
1165 Finternal.setZero();
1166 HE.setZero();
1167 GDEPSP.setZero();
1168 KALPHA.setZero();
1169
1170 // Enhanced Assumed Strain (EAS)
1171 T0.setZero();
1172 detJ0C = 0.0;
1173 T0DetJElementCenterForEAS(m_d0, T0, detJ0C);
1174 //== F_internal ==//
1175 // Choose constructors depending on m_isMooney
1176 Brick_ForceNumerical myformula =
1177 !m_isMooney ? Brick_ForceNumerical(&d, &m_d0, this, &T0, &detJ0C, &alpha_eas, &E, &v)
1178 : Brick_ForceNumerical(&d, &m_d0, this, &T0, &detJ0C, &alpha_eas);
1179 TempIntegratedResult.setZero();
1180 ChQuadrature::Integrate3D<ChVectorN<double, 330>>(
1181 TempIntegratedResult, // result of integration will go there
1182 myformula, // formula to integrate
1183 -1, // start of x
1184 1, // end of x
1185 -1, // start of y
1186 1, // end of y
1187 -1, // start of z
1188 1, // end of z
1189 2 // order of integration
1190 );
1191
1192 ///===============================================================//
1193 ///===TempIntegratedResult(0:23,1) -> InternalForce(24x1)=========//
1194 ///===TempIntegratedResult(24:32,1) -> HE(9x1) brick =========//
1195 ///===TempIntegratedResult(33:248,1) -> GDEPSP(9x24) =========//
1196 ///===TempIntegratedResult(249:329,1) -> KALPHA(9x9) =========//
1197 ///===============================================================//
1198 ChVectorN<double, 216> GDEPSPvec;
1199 ChVectorN<double, 81> KALPHAvec;
1200 Finternal = TempIntegratedResult.segment(0, 24);
1201 HE = TempIntegratedResult.segment(24, 9);
1202 GDEPSPvec = TempIntegratedResult.segment(33, 216);
1203 KALPHAvec = TempIntegratedResult.segment(249, 81);
1204 GDEPSP = Eigen::Map<ChMatrixNM<double, 9, 24>>(GDEPSPvec.data(), 9, 24);
1205 KALPHA = Eigen::Map<ChMatrixNM<double, 9, 9>>(KALPHAvec.data(), 9, 9);
1206 KALPHA1 = KALPHA;
1207
1208 if (m_flag_HE == NUMERICAL)
1209 break; // When numerical jacobian loop, no need to calculate HE
1210 count = count + 1;
1211 double norm_HE = HE.norm();
1212
1213 if (norm_HE < 0.00001) {
1214 fail = 0;
1215 } else {
1216 // Solve for ResidHE
1217 ResidHE = KALPHA1.colPivHouseholderQr().solve(HE);
1218 }
1219 if (m_flag_HE == ANALYTICAL && count > 2) {
1220 GetLog() << ie << " count " << count << " NormHE " << norm_HE << "\n";
1221 }
1222 }
1223 Fi = -Finternal;
1224 //== Stock_Alpha=================//
1225 if (m_flag_HE == ANALYTICAL) {
1226 SetStockAlpha(renewed_alpha_eas(0), renewed_alpha_eas(1), renewed_alpha_eas(2), renewed_alpha_eas(3),
1227 renewed_alpha_eas(4), renewed_alpha_eas(5), renewed_alpha_eas(6), renewed_alpha_eas(7),
1228 renewed_alpha_eas(8)); // this->
1229 }
1230 //== Jacobian Matrix for alpha ==//
1231 if (m_flag_HE == ANALYTICAL) {
1232 ChMatrixNM<double, 9, 9> INV_KALPHA;
1233 ChMatrixNM<double, 24, 24> stock_jac_EAS_elem;
1234
1235 for (int ii = 0; ii < 9; ii++) {
1236 ChVectorN<double, 9> DAMMY_vec;
1237 DAMMY_vec.setZero();
1238 DAMMY_vec(ii) = 1.0;
1239 INV_KALPHA.col(ii) = KALPHA.colPivHouseholderQr().solve(DAMMY_vec);
1240 }
1241 stock_jac_EAS_elem = GDEPSP.transpose() * INV_KALPHA * GDEPSP;
1242 m_stock_jac_EAS = stock_jac_EAS_elem;
1243 }
1244 } else {
1245 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1246
1247 ChVectorN<double, 906> TempIntegratedResult;
1248 ChVectorN<double, 24> Finternal;
1249 // Enhanced Assumed Strain (EAS)
1250 ChMatrixNM<double, 6, 6> T0;
1251 ChVectorN<double, 9> HE;
1252 ChMatrixNM<double, 9, 24> GDEPSP;
1253 ChMatrixNM<double, 9, 9> KALPHA;
1254 ChMatrixNM<double, 24, 24> KTE;
1255 ChMatrixNM<double, 9, 9> KALPHA1;
1256 ChVectorN<double, 9> ResidHE;
1257 double detJ0C;
1258 ChVectorN<double, 9> alpha_eas;
1259 ChVectorN<double, 9> renewed_alpha_eas;
1260 ChVectorN<double, 9> previous_alpha;
1261
1262 previous_alpha = m_stock_alpha_EAS;
1263 alpha_eas = previous_alpha;
1264 ResidHE.setZero();
1265 int count = 0;
1266 int fail = 1;
1267 // Loop to obtain convergence in EAS internal parameters alpha
1268 // This loops call ChQuadrature::Integrate3D on MyAnalyticalForce,
1269 // which calculates the Jacobian at every iteration of each time step
1270 int iteralpha = 0; // Counts number of iterations
1271 while (fail == 1) {
1272 iteralpha++;
1273 alpha_eas = alpha_eas - ResidHE;
1274 renewed_alpha_eas = alpha_eas;
1275
1276 Finternal.setZero(); // Internal force vector
1277 HE.setZero(); // Internal force vector from EAS
1278 GDEPSP.setZero(); // Jacobian of EAS forces w.r.t. coordinates
1279 KALPHA.setZero(); // Jacobian of EAS forces w.r.t. EAS internal parameters
1280
1281 // Enhanced Assumed Strain (EAS)
1282 T0.setZero();
1283 detJ0C = 0.0;
1284 T0DetJElementCenterForEAS(m_d0, T0, detJ0C);
1285
1286 //== F_internal ==//
1287 Brick_ForceAnalytical myformula =
1288 !m_isMooney ? Brick_ForceAnalytical(&d, &m_d0, this, &T0, &detJ0C, &alpha_eas, &E, &v)
1289 : Brick_ForceAnalytical(&d, &m_d0, this, &T0, &detJ0C, &alpha_eas);
1290 TempIntegratedResult.setZero();
1291 ChQuadrature::Integrate3D<ChVectorN<double, 906>>(
1292 TempIntegratedResult, // result of integration will go there
1293 myformula, // formula to integrate
1294 -1, // start of x
1295 1, // end of x
1296 -1, // start of y
1297 1, // end of y
1298 -1, // start of z
1299 1, // end of z
1300 2 // order of integration
1301 );
1302
1303 // ///===============================================================//
1304 // ///===TempIntegratedResult(0:23,1) -> InternalForce(24x1)=========//
1305 // ///===TempIntegratedResult(24:28,1) -> HE(5x1) =========//
1306 // ///===TempIntegratedResult(29:148,1) -> GDEPSP(5x24) =========//
1307 // ///===TempIntegratedResult(149:173,1) -> KALPHA(5x5) =========//
1308 // ///===TempIntegratedResult(174:749,1) -> Stiffness Matrix(24x24) =//
1309 // ///===============================================================//
1310 ChVectorN<double, 216> GDEPSPvec;
1311 ChVectorN<double, 81> KALPHAvec;
1312 ChVectorN<double, 576> JACvec;
1313 Finternal = TempIntegratedResult.segment(0, 24);
1314 HE = TempIntegratedResult.segment(24, 9);
1315 GDEPSPvec = TempIntegratedResult.segment(33, 216);
1316 KALPHAvec = TempIntegratedResult.segment(249, 81);
1317 JACvec = TempIntegratedResult.segment(330, 576);
1318 for (int i = 0; i < 9; i++) {
1319 for (int j = 0; j < 24; j++) {
1320 GDEPSP(i, j) = GDEPSPvec(i * 24 + j);
1321 }
1322 }
1323 // GDEPSP = GDEPSPvec;
1324 for (int i = 0; i < 9; i++) {
1325 for (int j = 0; j < 9; j++) {
1326 KALPHA(i, j) = KALPHAvec(i * 9 + j);
1327 }
1328 }
1329 // KALPHA = KALPHAvec;
1330 for (int i = 0; i < 24; i++) {
1331 for (int j = 0; j < 24; j++) {
1332 KTE(i, j) = JACvec(i * 24 + j);
1333 }
1334 }
1335 // KTE = JACvec;
1336
1337 // Calculation of the element Jacobian for implicit integrator
1338 // KTE and stock_jac_EAS_elem.
1339 KALPHA1 = KALPHA;
1340 if (m_flag_HE == NUMERICAL)
1341 break; // When numerical jacobian loop, no need to calculate HE
1342 count = count + 1;
1343 double norm_HE = HE.norm();
1344 if (norm_HE < 0.00001) {
1345 fail = 0;
1346 } else {
1347 // Solve for ResidHE
1348 ResidHE = KALPHA1.colPivHouseholderQr().solve(HE);
1349 }
1350 } // end of while
1351 Fi = -Finternal;
1352 ////== Stock_Alpha=================//
1353 if (m_flag_HE == ANALYTICAL) {
1354 SetStockAlpha(renewed_alpha_eas(0), renewed_alpha_eas(1), renewed_alpha_eas(2), renewed_alpha_eas(3),
1355 renewed_alpha_eas(4), renewed_alpha_eas(5), renewed_alpha_eas(6), renewed_alpha_eas(7),
1356 renewed_alpha_eas(8)); // this->
1357 }
1358 ////== Jacobian Matrix for alpha ==//
1359 if (m_flag_HE == ANALYTICAL) {
1360 ChMatrixNM<double, 9, 9> INV_KALPHA;
1361 ChMatrixNM<double, 24, 24> stock_jac_EAS_elem;
1362
1363 for (int ii = 0; ii < 9; ii++) {
1364 ChVectorN<double, 9> DAMMY_vec;
1365 DAMMY_vec.setZero();
1366 DAMMY_vec(ii) = 1.0;
1367 INV_KALPHA.col(ii) = KALPHA.colPivHouseholderQr().solve(DAMMY_vec);
1368 }
1369 stock_jac_EAS_elem = GDEPSP.transpose() * INV_KALPHA * GDEPSP;
1370 m_stock_KTE = KTE;
1371 m_stock_jac_EAS = stock_jac_EAS_elem;
1372 }
1373 } // end of else for numerical or analytical
1374 }
1375
1376 // -----------------------------------------------------------------------------
1377
ShapeFunctions(ShapeVector & N,double x,double y,double z)1378 void ChElementHexaANCF_3813::ShapeFunctions(ShapeVector& N, double x, double y, double z) {
1379 N(0) = 0.125 * (1.0 - x) * (1.0 - y) * (1.0 - z);
1380 N(1) = 0.125 * (1.0 + x) * (1.0 - y) * (1.0 - z);
1381 N(2) = 0.125 * (1.0 + x) * (1.0 + y) * (1.0 - z);
1382 N(3) = 0.125 * (1.0 - x) * (1.0 + y) * (1.0 - z);
1383 N(4) = 0.125 * (1.0 - x) * (1.0 - y) * (1.0 + z);
1384 N(5) = 0.125 * (1.0 + x) * (1.0 - y) * (1.0 + z);
1385 N(6) = 0.125 * (1.0 + x) * (1.0 + y) * (1.0 + z);
1386 N(7) = 0.125 * (1.0 - x) * (1.0 + y) * (1.0 + z);
1387 }
1388
1389 // -----------------------------------------------------------------------------
1390
ShapeFunctionsDerivativeX(ShapeVector & Nx,double x,double y,double z)1391 void ChElementHexaANCF_3813::ShapeFunctionsDerivativeX(ShapeVector& Nx, double x, double y, double z) {
1392 double a = GetLengthX();
1393
1394 Nx(0) = 2.0 / a * 0.125 * (-1.0) * (1.0 - y) * (1.0 - z);
1395 Nx(1) = 2.0 / a * 0.125 * (1.0) * (1.0 - y) * (1.0 - z);
1396 Nx(2) = 2.0 / a * 0.125 * (1.0) * (1.0 + y) * (1.0 - z);
1397 Nx(3) = 2.0 / a * 0.125 * (-1.0) * (1.0 + y) * (1.0 - z);
1398 Nx(4) = 2.0 / a * 0.125 * (-1.0) * (1.0 - y) * (1.0 + z);
1399 Nx(5) = 2.0 / a * 0.125 * (1.0) * (1.0 - y) * (1.0 + z);
1400 Nx(6) = 2.0 / a * 0.125 * (1.0) * (1.0 + y) * (1.0 + z);
1401 Nx(7) = 2.0 / a * 0.125 * (-1.0) * (1.0 + y) * (1.0 + z);
1402 }
1403
1404 // -----------------------------------------------------------------------------
1405
ShapeFunctionsDerivativeY(ShapeVector & Ny,double x,double y,double z)1406 void ChElementHexaANCF_3813::ShapeFunctionsDerivativeY(ShapeVector& Ny, double x, double y, double z) {
1407 double b = GetLengthY();
1408
1409 Ny(0) = 2.0 / b * 0.125 * (1.0 - x) * (-1.0) * (1.0 - z);
1410 Ny(1) = 2.0 / b * 0.125 * (1.0 + x) * (-1.0) * (1.0 - z);
1411 Ny(2) = 2.0 / b * 0.125 * (1.0 + x) * (1.0) * (1.0 - z);
1412 Ny(3) = 2.0 / b * 0.125 * (1.0 - x) * (1.0) * (1.0 - z);
1413 Ny(4) = 2.0 / b * 0.125 * (1.0 - x) * (-1.0) * (1.0 + z);
1414 Ny(5) = 2.0 / b * 0.125 * (1.0 + x) * (-1.0) * (1.0 + z);
1415 Ny(6) = 2.0 / b * 0.125 * (1.0 + x) * (1.0) * (1.0 + z);
1416 Ny(7) = 2.0 / b * 0.125 * (1.0 - x) * (1.0) * (1.0 + z);
1417 }
1418
1419 // -----------------------------------------------------------------------------
1420
ShapeFunctionsDerivativeZ(ShapeVector & Nz,double x,double y,double z)1421 void ChElementHexaANCF_3813::ShapeFunctionsDerivativeZ(ShapeVector& Nz, double x, double y, double z) {
1422 double c = GetLengthZ();
1423
1424 Nz(0) = 2.0 / c * 0.125 * (1.0 - x) * (1.0 - y) * (-1.0);
1425 Nz(1) = 2.0 / c * 0.125 * (1.0 + x) * (1.0 - y) * (-1.0);
1426 Nz(2) = 2.0 / c * 0.125 * (1.0 + x) * (1.0 + y) * (-1.0);
1427 Nz(3) = 2.0 / c * 0.125 * (1.0 - x) * (1.0 + y) * (-1.0);
1428 Nz(4) = 2.0 / c * 0.125 * (1.0 - x) * (1.0 - y) * (1.0);
1429 Nz(5) = 2.0 / c * 0.125 * (1.0 + x) * (1.0 - y) * (1.0);
1430 Nz(6) = 2.0 / c * 0.125 * (1.0 + x) * (1.0 + y) * (1.0);
1431 Nz(7) = 2.0 / c * 0.125 * (1.0 - x) * (1.0 + y) * (1.0);
1432 }
1433
1434 // ----------------------------------------------------------------------------
1435
Update()1436 void ChElementHexaANCF_3813::Update() {
1437 // parent class update:
1438 ChElementGeneric::Update();
1439 }
1440
1441 // -----------------------------------------------------------------------------
1442
GetStateBlock(ChVectorDynamic<> & mD)1443 void ChElementHexaANCF_3813::GetStateBlock(ChVectorDynamic<>& mD) {
1444 mD.segment(0, 3) = m_nodes[0]->GetPos().eigen();
1445 mD.segment(3, 3) = m_nodes[1]->GetPos().eigen();
1446 mD.segment(6, 3) = m_nodes[2]->GetPos().eigen();
1447 mD.segment(9, 3) = m_nodes[3]->GetPos().eigen();
1448 mD.segment(12, 3) = m_nodes[4]->GetPos().eigen();
1449 mD.segment(15, 3) = m_nodes[5]->GetPos().eigen();
1450 mD.segment(18, 3) = m_nodes[6]->GetPos().eigen();
1451 mD.segment(21, 3) = m_nodes[7]->GetPos().eigen();
1452 }
1453
1454 // -----------------------------------------------------------------------------
1455
ComputeStiffnessMatrix()1456 void ChElementHexaANCF_3813::ComputeStiffnessMatrix() {
1457 bool use_numerical_differentiation = false;
1458
1459 if (use_numerical_differentiation) {
1460 double diff = 1e-8;
1461 ChVectorDynamic<> F0(24);
1462 ChVectorDynamic<> F1(24);
1463 ComputeInternalForces(F0);
1464 for (int inode = 0; inode < 8; ++inode) {
1465 m_nodes[inode]->pos.x() += diff;
1466 ComputeInternalForces(F1); // Flag=1 > Jacobian of internal force calculation
1467 m_StiffnessMatrix.col(0 + inode * 3) = (F0 - F1) * (1.0 / diff);
1468 m_nodes[inode]->pos.x() -= diff;
1469
1470 m_nodes[inode]->pos.y() += diff;
1471 ComputeInternalForces(F1);
1472 m_StiffnessMatrix.col(1 + inode * 3) = (F0 - F1) * (1.0 / diff);
1473 m_nodes[inode]->pos.y() -= diff;
1474
1475 m_nodes[inode]->pos.z() += diff;
1476 ComputeInternalForces(F1);
1477 m_StiffnessMatrix.col(2 + inode * 3) = (F0 - F1) * (1.0 / diff);
1478 m_nodes[inode]->pos.z() -= diff;
1479 }
1480 // flag_HE=0 is default
1481 m_StiffnessMatrix -= m_stock_jac_EAS; // For Enhanced Assumed Strain
1482 } else {
1483 // Put in m_StiffnessMatrix the values for the Jacobian already calculated in the computation of internal forces
1484 // Note that m_stock_KTE and m_stock_jac_EAS are updated at each iteration of the time step
1485 m_StiffnessMatrix = m_stock_KTE;
1486 m_StiffnessMatrix -= m_stock_jac_EAS;
1487 }
1488 }
1489
1490 // -----------------------------------------------------------------------------
1491
1492 class Brick_Mass : public ChIntegrable3D<ChMatrixNM<double, 24, 24>> {
1493 public:
1494 Brick_Mass(ChMatrixNM<double, 8, 3>* d0_, ChElementHexaANCF_3813* element_);
~Brick_Mass()1495 ~Brick_Mass() {}
1496
1497 private:
1498 ChElementHexaANCF_3813* element;
1499 ChMatrixNM<double, 8, 3>* d0; ///< Pointer to a matrix containing the element initial coordinates
1500 ChMatrixNM<double, 3, 24> S; ///< Sparse shape function matrix
1501 ChElementHexaANCF_3813::ShapeVector N; ///< Dense shape function vector
1502 ChElementHexaANCF_3813::ShapeVector Nx; ///< Dense shape function vector, X derivative
1503 ChElementHexaANCF_3813::ShapeVector Ny; ///< Dense shape function vector, Y derivative
1504 ChElementHexaANCF_3813::ShapeVector Nz; ///< Dense shape function vector, Z derivative
1505
1506 /// Evaluate the S'*S at point x
1507 virtual void Evaluate(ChMatrixNM<double, 24, 24>& result, const double x, const double y, const double z) override;
1508 };
1509
Brick_Mass(ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_)1510 Brick_Mass::Brick_Mass(ChMatrixNM<double, 8, 3>* d0_, ChElementHexaANCF_3813* element_) : element(element_), d0(d0_) {
1511 S.setZero();
1512 }
1513
Evaluate(ChMatrixNM<double,24,24> & result,const double x,const double y,const double z)1514 void Brick_Mass::Evaluate(ChMatrixNM<double, 24, 24>& result, const double x, const double y, const double z) {
1515 element->ShapeFunctions(N, x, y, z);
1516 element->ShapeFunctionsDerivativeX(Nx, x, y, z);
1517 element->ShapeFunctionsDerivativeY(Ny, x, y, z);
1518 element->ShapeFunctionsDerivativeZ(Nz, x, y, z);
1519
1520 // S=[N1*eye(3) N2*eye(3) N3*eye(3) N4*eye(3)...]
1521 for (int i = 0; i < 8; i++) {
1522 S(0, 3 * i + 0) = N(i);
1523 S(1, 3 * i + 1) = N(i);
1524 S(2, 3 * i + 2) = N(i);
1525 }
1526
1527 ChMatrixNM<double, 3, 3> rd0;
1528 rd0.col(0) = (*d0).transpose() * Nx.transpose();
1529 rd0.col(1) = (*d0).transpose() * Ny.transpose();
1530 rd0.col(2) = (*d0).transpose() * Nz.transpose();
1531 double detJ0 = rd0.determinant();
1532
1533 // Perform r = S'*S
1534 double factor = detJ0 * (element->GetLengthX() / 2) * (element->GetLengthY() / 2) * (element->GetLengthZ() / 2);
1535 result = factor * S.transpose() * S;
1536 }
1537
ComputeMassMatrix()1538 void ChElementHexaANCF_3813::ComputeMassMatrix() {
1539 double rho = m_Material->Get_density();
1540 Brick_Mass myformula(&m_d0, this);
1541 m_MassMatrix.setZero();
1542 ChQuadrature::Integrate3D<ChMatrixNM<double, 24, 24>>(m_MassMatrix, // result of integration will go there
1543 myformula, // formula to integrate
1544 -1, // start of x
1545 1, // end of x
1546 -1, // start of y
1547 1, // end of y
1548 -1, // start of z
1549 1, // end of z
1550 2 // order of integration
1551 );
1552
1553 m_MassMatrix *= rho;
1554 }
1555
1556 // -----------------------------------------------------------------------------
1557
1558 // Class to calculate the gravity forces of a brick element
1559 class BrickGravity : public ChIntegrable3D<ChVectorN<double, 8>> {
1560 public:
1561 BrickGravity(ChMatrixNM<double, 8, 3>* d0_, ChElementHexaANCF_3813* element_);
~BrickGravity()1562 ~BrickGravity() {}
1563
1564 private:
1565 ChElementHexaANCF_3813* element;
1566 ChMatrixNM<double, 8, 3>* d0; // Pointer to a matrix containing the element initial coordinates
1567 ChMatrixNM<double, 3, 24> S; // Sparse shape function matrix
1568 ChElementHexaANCF_3813::ShapeVector N; // Dense shape function vector
1569 ChElementHexaANCF_3813::ShapeVector Nx; // Dense shape function vector, X derivative
1570 ChElementHexaANCF_3813::ShapeVector Ny; // Dense shape function vector, Y derivative
1571 ChElementHexaANCF_3813::ShapeVector Nz; // Dense shape function vector, Z derivative
1572
1573 virtual void Evaluate(ChVectorN<double, 8>& result, const double x, const double y, const double z) override;
1574 };
1575
BrickGravity(ChMatrixNM<double,8,3> * d0_,ChElementHexaANCF_3813 * element_)1576 BrickGravity::BrickGravity(ChMatrixNM<double, 8, 3>* d0_, ChElementHexaANCF_3813* element_)
1577 : element(element_), d0(d0_) {}
1578
Evaluate(ChVectorN<double,8> & result,const double x,const double y,const double z)1579 void BrickGravity::Evaluate(ChVectorN<double, 8>& result, const double x, const double y, const double z) {
1580 element->ShapeFunctions(N, x, y, z);
1581 element->ShapeFunctionsDerivativeX(Nx, x, y, z);
1582 element->ShapeFunctionsDerivativeY(Ny, x, y, z);
1583 element->ShapeFunctionsDerivativeZ(Nz, x, y, z);
1584
1585 // Weights for Gaussian integration
1586 double wx2 = (element->GetLengthX()) / 2;
1587 double wy2 = (element->GetLengthY()) / 2;
1588 double wz2 = (element->GetLengthZ()) / 2;
1589
1590 ChMatrixNM<double, 3, 3> rd0;
1591 rd0.col(0) = (*d0).transpose() * Nx.transpose();
1592 rd0.col(1) = (*d0).transpose() * Ny.transpose();
1593 rd0.col(2) = (*d0).transpose() * Nz.transpose();
1594 double detJ0 = rd0.determinant();
1595
1596 result = detJ0 * wx2 * wy2 * wz2 * N.transpose();
1597 }
1598
ComputeGravityForceScale()1599 void ChElementHexaANCF_3813::ComputeGravityForceScale() {
1600 BrickGravity myformula1(&m_d0, this);
1601 m_GravForceScale.setZero();
1602 ChQuadrature::Integrate3D<ChVectorN<double, 8>>(m_GravForceScale, // result of integration will go there
1603 myformula1, // formula to integrate
1604 -1, 1, // limits in x direction
1605 -1, 1, // limits in y direction
1606 -1, 1, // limits in z direction
1607 2 // order of integration
1608 );
1609
1610 m_GravForceScale *= m_Material->Get_density();
1611 }
1612
1613 // Compute the generalized force vector due to gravity
ComputeGravityForces(ChVectorDynamic<> & Fg,const ChVector<> & G_acc)1614 void ChElementHexaANCF_3813::ComputeGravityForces(ChVectorDynamic<>& Fg, const ChVector<>& G_acc) {
1615 assert(Fg.size() == 24);
1616
1617 // Calculate and add the generalized force due to gravity to the generalized internal force vector for the element.
1618 // The generalized force due to gravity could be computed once prior to the start of the simulation if gravity was
1619 // assumed constant throughout the entire simulation. However, this implementation assumes that the acceleration
1620 // due to gravity, while a constant for the entire system, can change from step to step which could be useful for
1621 // gravity loaded units tests as an example. The generalized force due to gravity is calculated in compact matrix
1622 // form and is pre-mapped to the desired vector format
1623 Eigen::Map<ChMatrixNM<double, 8, 3>> GravForceCompact(Fg.data(), 8, 3);
1624 GravForceCompact = m_GravForceScale * G_acc.eigen().transpose();
1625 }
1626
1627 // -----------------------------------------------------------------------------
1628
SetupInitial(ChSystem * system)1629 void ChElementHexaANCF_3813::SetupInitial(ChSystem* system) {
1630 // Compute gravitational forces
1631 ComputeGravityForceScale();
1632 // Compute mass matrix
1633 ComputeMassMatrix();
1634 // initial EAS parameters
1635 m_stock_jac_EAS.setZero();
1636 // Compute stiffness matrix
1637 // (this is not constant in ANCF and will be called automatically many times by ComputeKRMmatricesGlobal()
1638 // when the solver will run, yet maybe nice to privide an initial nonzero value)
1639 ComputeStiffnessMatrix();
1640 }
1641
1642 // -----------------------------------------------------------------------------
1643
ComputeKRMmatricesGlobal(ChMatrixRef H,double Kfactor,double Rfactor,double Mfactor)1644 void ChElementHexaANCF_3813::ComputeKRMmatricesGlobal(ChMatrixRef H, double Kfactor, double Rfactor, double Mfactor) {
1645 assert((H.rows() == 24) && (H.cols() == 24));
1646
1647 // Compute global stiffness matrix:
1648 ComputeStiffnessMatrix();
1649
1650 // For K stiffness matrix and R matrix: scale by factors
1651 // because [R] = r*[K] , so kf*[K]+rf*[R] = (kf+rf*r)*[K]
1652 double kr_factor = Kfactor + Rfactor * m_Material->Get_RayleighDampingK();
1653
1654 // Paste scaled K stiffness matrix and R matrix in resulting H and add scaled mass matrix.
1655 H.block(0, 0, 24, 24) = kr_factor * m_StiffnessMatrix + Mfactor * m_MassMatrix;
1656 }
1657
1658 // -----------------------------------------------------------------------------
1659
T0DetJElementCenterForEAS(ChMatrixNM<double,8,3> & d0,ChMatrixNM<double,6,6> & T0,double & detJ0C)1660 void ChElementHexaANCF_3813::T0DetJElementCenterForEAS(ChMatrixNM<double, 8, 3>& d0,
1661 ChMatrixNM<double, 6, 6>& T0,
1662 double& detJ0C) {
1663 ShapeVector Nx;
1664 ShapeVector Ny;
1665 ShapeVector Nz;
1666 ShapeFunctionsDerivativeX(Nx, 0, 0, 0);
1667 ShapeFunctionsDerivativeY(Ny, 0, 0, 0);
1668 ShapeFunctionsDerivativeZ(Nz, 0, 0, 0);
1669
1670 ChMatrixNM<double, 3, 3> rd0;
1671 rd0.col(0) = d0.transpose() * Nx.transpose();
1672 rd0.col(1) = d0.transpose() * Ny.transpose();
1673 rd0.col(2) = d0.transpose() * Nz.transpose();
1674 detJ0C = rd0.determinant();
1675
1676 // Transformation : Orthogonal transformation (A and J) ////
1677 ChVector<double> G1;
1678 ChVector<double> G2;
1679 ChVector<double> G3;
1680 ChVector<double> G1xG2;
1681 G1[0] = rd0(0, 0);
1682 G2[0] = rd0(0, 1);
1683 G3[0] = rd0(0, 2);
1684 G1[1] = rd0(1, 0);
1685 G2[1] = rd0(1, 1);
1686 G3[1] = rd0(1, 2);
1687 G1[2] = rd0(2, 0);
1688 G2[2] = rd0(2, 1);
1689 G3[2] = rd0(2, 2);
1690 G1xG2.Cross(G1, G2);
1691
1692 // Tangent Frame
1693 ChVector<> A1 = G1 / sqrt(G1[0] * G1[0] + G1[1] * G1[1] + G1[2] * G1[2]);
1694 ChVector<> A3 = G1xG2 / sqrt(G1xG2[0] * G1xG2[0] + G1xG2[1] * G1xG2[1] + G1xG2[2] * G1xG2[2]);
1695 ChVector<> A2 = A3.Cross(A1);
1696 double theta = 0.0;
1697 ChVector<> AA1 = A1 * cos(theta) + A2 * sin(theta);
1698 ChVector<> AA2 = -A1 * sin(theta) + A2 * cos(theta);
1699 ChVector<> AA3 = A3;
1700
1701 // Beta
1702 ChMatrixNM<double, 3, 3> j0 = rd0.inverse();
1703 ChVector<double> j01;
1704 ChVector<double> j02;
1705 ChVector<double> j03;
1706 j01[0] = j0(0, 0);
1707 j02[0] = j0(1, 0);
1708 j03[0] = j0(2, 0);
1709 j01[1] = j0(0, 1);
1710 j02[1] = j0(1, 1);
1711 j03[1] = j0(2, 1);
1712 j01[2] = j0(0, 2);
1713 j02[2] = j0(1, 2);
1714 j03[2] = j0(2, 2);
1715 ChVectorN<double, 9> beta;
1716 beta(0) = Vdot(AA1, j01);
1717 beta(1) = Vdot(AA2, j01);
1718 beta(2) = Vdot(AA3, j01);
1719 beta(3) = Vdot(AA1, j02);
1720 beta(4) = Vdot(AA2, j02);
1721 beta(5) = Vdot(AA3, j02);
1722 beta(6) = Vdot(AA1, j03);
1723 beta(7) = Vdot(AA2, j03);
1724 beta(8) = Vdot(AA3, j03);
1725
1726 T0(0, 0) = pow(beta(0), 2);
1727 T0(1, 0) = pow(beta(1), 2);
1728 T0(2, 0) = 2.0 * beta(0) * beta(1);
1729 T0(3, 0) = pow(beta(2), 2);
1730 T0(4, 0) = 2.0 * beta(0) * beta(2);
1731 T0(5, 0) = 2.0 * beta(1) * beta(2);
1732
1733 T0(0, 1) = pow(beta(3), 2);
1734 T0(1, 1) = pow(beta(4), 2);
1735 T0(2, 1) = 2.0 * beta(3) * beta(4);
1736 T0(3, 1) = pow(beta(5), 2);
1737 T0(4, 1) = 2.0 * beta(3) * beta(5);
1738 T0(5, 1) = 2.0 * beta(4) * beta(5);
1739
1740 T0(0, 2) = beta(0) * beta(3);
1741 T0(1, 2) = beta(1) * beta(4);
1742 T0(2, 2) = beta(0) * beta(4) + beta(1) * beta(3);
1743 T0(3, 2) = beta(2) * beta(5);
1744 T0(4, 2) = beta(0) * beta(5) + beta(2) * beta(3);
1745 T0(5, 2) = beta(2) * beta(4) + beta(1) * beta(5);
1746
1747 T0(0, 3) = pow(beta(6), 2);
1748 T0(1, 3) = pow(beta(7), 2);
1749 T0(2, 3) = 2.0 * beta(6) * beta(7);
1750 T0(3, 3) = pow(beta(8), 2);
1751 T0(4, 3) = 2.0 * beta(6) * beta(8);
1752 T0(5, 3) = 2.0 * beta(7) * beta(8);
1753
1754 T0(0, 4) = beta(0) * beta(6);
1755 T0(1, 4) = beta(1) * beta(7);
1756 T0(2, 4) = beta(0) * beta(7) + beta(6) * beta(1);
1757 T0(3, 4) = beta(2) * beta(8);
1758 T0(4, 4) = beta(0) * beta(8) + beta(2) * beta(6);
1759 T0(5, 4) = beta(1) * beta(8) + beta(2) * beta(7);
1760
1761 T0(0, 5) = beta(3) * beta(6);
1762 T0(1, 5) = beta(4) * beta(7);
1763 T0(2, 5) = beta(3) * beta(7) + beta(4) * beta(6);
1764 T0(3, 5) = beta(5) * beta(8);
1765 T0(4, 5) = beta(3) * beta(8) + beta(6) * beta(5);
1766 T0(5, 5) = beta(4) * beta(8) + beta(5) * beta(7);
1767 }
1768
Basis_M(ChMatrixNM<double,6,9> & M,double x,double y,double z)1769 void ChElementHexaANCF_3813::Basis_M(ChMatrixNM<double, 6, 9>& M, double x, double y, double z) {
1770 M.setZero();
1771 M(0, 0) = x;
1772 M(1, 1) = y;
1773 M(2, 2) = x;
1774 M(2, 3) = y;
1775 M(3, 4) = z;
1776 M(4, 5) = x;
1777 M(4, 6) = z;
1778 M(5, 7) = y;
1779 M(5, 8) = z;
1780 }
1781
1782 // -----------------------------------------------------------------------------
1783
LoadableGetStateBlock_x(int block_offset,ChState & mD)1784 void ChElementHexaANCF_3813::LoadableGetStateBlock_x(int block_offset, ChState& mD) {
1785 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos().eigen();
1786 mD.segment(block_offset + 3, 3) = m_nodes[1]->GetPos().eigen();
1787 mD.segment(block_offset + 6, 3) = m_nodes[2]->GetPos().eigen();
1788 mD.segment(block_offset + 9, 3) = m_nodes[3]->GetPos().eigen();
1789 mD.segment(block_offset + 12, 3) = m_nodes[4]->GetPos().eigen();
1790 mD.segment(block_offset + 15, 3) = m_nodes[5]->GetPos().eigen();
1791 mD.segment(block_offset + 18, 3) = m_nodes[6]->GetPos().eigen();
1792 mD.segment(block_offset + 21, 3) = m_nodes[7]->GetPos().eigen();
1793 }
1794
LoadableGetStateBlock_w(int block_offset,ChStateDelta & mD)1795 void ChElementHexaANCF_3813::LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) {
1796 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos_dt().eigen();
1797 mD.segment(block_offset + 3, 3) = m_nodes[1]->GetPos_dt().eigen();
1798 mD.segment(block_offset + 6, 3) = m_nodes[2]->GetPos_dt().eigen();
1799 mD.segment(block_offset + 9, 3) = m_nodes[3]->GetPos_dt().eigen();
1800 mD.segment(block_offset + 12, 3) = m_nodes[4]->GetPos_dt().eigen();
1801 mD.segment(block_offset + 15, 3) = m_nodes[5]->GetPos_dt().eigen();
1802 mD.segment(block_offset + 18, 3) = m_nodes[6]->GetPos_dt().eigen();
1803 mD.segment(block_offset + 21, 3) = m_nodes[7]->GetPos_dt().eigen();
1804 }
1805
LoadableStateIncrement(const unsigned int off_x,ChState & x_new,const ChState & x,const unsigned int off_v,const ChStateDelta & Dv)1806 void ChElementHexaANCF_3813::LoadableStateIncrement(const unsigned int off_x,
1807 ChState& x_new,
1808 const ChState& x,
1809 const unsigned int off_v,
1810 const ChStateDelta& Dv) {
1811 for (int i = 0; i < 8; ++i) {
1812 this->m_nodes[i]->NodeIntStateIncrement(off_x + 3 * 1, x_new, x, off_v + 3 * i, Dv);
1813 }
1814 }
1815
LoadableGetVariables(std::vector<ChVariables * > & mvars)1816 void ChElementHexaANCF_3813::LoadableGetVariables(std::vector<ChVariables*>& mvars) {
1817 for (int i = 0; i < m_nodes.size(); ++i)
1818 mvars.push_back(&this->m_nodes[i]->Variables());
1819 }
1820
1821 // -----------------------------------------------------------------------------
1822
ComputeNF(const double U,const double V,const double W,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1823 void ChElementHexaANCF_3813::ComputeNF(
1824 const double U, // parametric coordinate in volume
1825 const double V, // parametric coordinate in volume
1826 const double W, // parametric coordinate in volume
1827 ChVectorDynamic<>& Qi, // Return result of N'*F here, maybe with offset block_offset
1828 double& detJ, // Return det[J] here
1829 const ChVectorDynamic<>& F, // Input F vector, size is = n.field coords.
1830 ChVectorDynamic<>* state_x, // if != 0, update state (pos. part) to this, then evaluate Q
1831 ChVectorDynamic<>* state_w // if != 0, update state (speed part) to this, then evaluate Q
1832 ) {
1833 // this->ComputeNF(U, V, Qi, detJ, F, state_x, state_w);
1834 ShapeVector N;
1835 ShapeVector Nx;
1836 ShapeVector Ny;
1837 ShapeVector Nz;
1838 ShapeFunctions(N, U, V, W); // evaluate shape functions (in compressed vector)
1839 ShapeFunctionsDerivativeX(Nx, U, V, W);
1840 ShapeFunctionsDerivativeY(Ny, U, V, W);
1841 ShapeFunctionsDerivativeZ(Nz, U, V, W);
1842
1843 ChMatrixNM<double, 3, 3> rd0;
1844 rd0.col(0) = m_d0.transpose() * Nx.transpose();
1845 rd0.col(1) = m_d0.transpose() * Ny.transpose();
1846 rd0.col(2) = m_d0.transpose() * Nz.transpose();
1847 detJ = rd0.determinant();
1848 detJ *= this->GetLengthX() * this->GetLengthY() * this->GetLengthZ() / 8.0;
1849
1850 for (int i = 0; i < 8; i++) {
1851 Qi.segment(3 * i, 3) = N(i) * F.segment(0, 3);
1852 }
1853 }
1854
1855 } // end namespace fea
1856 } // end namespace chrono
1857