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, Milad Rakhsha, Antonio Recuero, Radu Serban
13 // =============================================================================
14 // ANCF laminated shell element with four nodes.
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 //// - remove unecessary initializations to zero
22
23 #include <cmath>
24
25 #include "chrono/fea/ChElementShellANCF_3423.h"
26 #include "chrono/core/ChException.h"
27 #include "chrono/core/ChQuadrature.h"
28 #include "chrono/physics/ChSystem.h"
29
30 namespace chrono {
31 namespace fea {
32
33 // ------------------------------------------------------------------------------
34 // Static variables
35 // ------------------------------------------------------------------------------
36 const double ChElementShellANCF_3423::m_toleranceEAS = 1e-5;
37 const int ChElementShellANCF_3423::m_maxIterationsEAS = 100;
38
39 // ------------------------------------------------------------------------------
40 // Constructor
41 // ------------------------------------------------------------------------------
42
ChElementShellANCF_3423()43 ChElementShellANCF_3423::ChElementShellANCF_3423() : m_numLayers(0), m_lenX(0), m_lenY(0), m_thickness(0), m_Alpha(0) {
44 m_nodes.resize(4);
45 }
46
47 // ------------------------------------------------------------------------------
48 // Set element nodes
49 // ------------------------------------------------------------------------------
50
SetNodes(std::shared_ptr<ChNodeFEAxyzD> nodeA,std::shared_ptr<ChNodeFEAxyzD> nodeB,std::shared_ptr<ChNodeFEAxyzD> nodeC,std::shared_ptr<ChNodeFEAxyzD> nodeD)51 void ChElementShellANCF_3423::SetNodes(std::shared_ptr<ChNodeFEAxyzD> nodeA,
52 std::shared_ptr<ChNodeFEAxyzD> nodeB,
53 std::shared_ptr<ChNodeFEAxyzD> nodeC,
54 std::shared_ptr<ChNodeFEAxyzD> nodeD) {
55 assert(nodeA);
56 assert(nodeB);
57 assert(nodeC);
58 assert(nodeD);
59
60 m_nodes[0] = nodeA;
61 m_nodes[1] = nodeB;
62 m_nodes[2] = nodeC;
63 m_nodes[3] = nodeD;
64 std::vector<ChVariables*> mvars;
65 mvars.push_back(&m_nodes[0]->Variables());
66 mvars.push_back(&m_nodes[0]->Variables_D());
67 mvars.push_back(&m_nodes[1]->Variables());
68 mvars.push_back(&m_nodes[1]->Variables_D());
69 mvars.push_back(&m_nodes[2]->Variables());
70 mvars.push_back(&m_nodes[2]->Variables_D());
71 mvars.push_back(&m_nodes[3]->Variables());
72 mvars.push_back(&m_nodes[3]->Variables_D());
73 Kmatr.SetVariables(mvars);
74
75 // Initial positions and slopes of the element nodes
76 CalcCoordMatrix(m_d0);
77 m_d0d0T = m_d0 * m_d0.transpose();
78 }
79
80 // -----------------------------------------------------------------------------
81 // Add a layer.
82 // -----------------------------------------------------------------------------
83
AddLayer(double thickness,double theta,std::shared_ptr<ChMaterialShellANCF> material)84 void ChElementShellANCF_3423::AddLayer(double thickness, double theta, std::shared_ptr<ChMaterialShellANCF> material) {
85 m_layers.push_back(Layer(this, thickness, theta, material));
86 }
87
88 // -----------------------------------------------------------------------------
89 // Interface to ChElementBase base class
90 // -----------------------------------------------------------------------------
91
92 // Initial element setup.
SetupInitial(ChSystem * system)93 void ChElementShellANCF_3423::SetupInitial(ChSystem* system) {
94 // Perform layer initialization and accumulate element thickness.
95 m_numLayers = m_layers.size();
96 m_thickness = 0;
97 for (size_t kl = 0; kl < m_numLayers; kl++) {
98 m_layers[kl].SetupInitial();
99 m_thickness += m_layers[kl].Get_thickness();
100 }
101
102 // Loop again over the layers and calculate the range for Gauss integration in the
103 // z direction (values in [-1,1]).
104 m_GaussZ.push_back(-1);
105 double z = 0;
106 for (size_t kl = 0; kl < m_numLayers; kl++) {
107 z += m_layers[kl].Get_thickness();
108 m_GaussZ.push_back(2 * z / m_thickness - 1);
109 }
110
111 // Reserve space for the EAS parameters and Jacobians.
112 m_alphaEAS.resize(m_numLayers);
113 m_KalphaEAS.resize(m_numLayers);
114 for (size_t i = 0; i < m_numLayers; i++) {
115 m_alphaEAS[i].setZero();
116 m_KalphaEAS[i].setZero();
117 }
118
119 // Cache the scaling factor (due to change of integration intervals)
120 m_GaussScaling = (m_lenX * m_lenY * m_thickness) / 8;
121
122 // Compute mass matrix and gravitational force scaling term (constant)
123 ComputeMassMatrix();
124 ComputeGravityForceScale();
125 }
126
127 // State update.
Update()128 void ChElementShellANCF_3423::Update() {
129 ChElementGeneric::Update();
130 }
131
132 // Fill the D vector with the current field values at the element nodes.
GetStateBlock(ChVectorDynamic<> & mD)133 void ChElementShellANCF_3423::GetStateBlock(ChVectorDynamic<>& mD) {
134 mD.segment(0, 3) = m_nodes[0]->GetPos().eigen();
135 mD.segment(3, 3) = m_nodes[0]->GetD().eigen();
136 mD.segment(6, 3) = m_nodes[1]->GetPos().eigen();
137 mD.segment(9, 3) = m_nodes[1]->GetD().eigen();
138 mD.segment(12, 3) = m_nodes[2]->GetPos().eigen();
139 mD.segment(15, 3) = m_nodes[2]->GetD().eigen();
140 mD.segment(18, 3) = m_nodes[3]->GetPos().eigen();
141 mD.segment(21, 3) = m_nodes[3]->GetD().eigen();
142 }
143
144 // Calculate the global matrix H as a linear combination of K, R, and M:
145 // H = Mfactor * [M] + Kfactor * [K] + Rfactor * [R]
ComputeKRMmatricesGlobal(ChMatrixRef H,double Kfactor,double Rfactor,double Mfactor)146 void ChElementShellANCF_3423::ComputeKRMmatricesGlobal(ChMatrixRef H, double Kfactor, double Rfactor, double Mfactor) {
147 assert((H.rows() == 24) && (H.cols() == 24));
148
149 // Calculate the linear combination Kfactor*[K] + Rfactor*[R]
150 ComputeInternalJacobians(Kfactor, Rfactor);
151
152 // Load Jac + Mfactor*[M] into H
153 for (int i = 0; i < 24; i++)
154 for (int j = 0; j < 24; j++)
155 H(i, j) = m_JacobianMatrix(i, j) + Mfactor * m_MassMatrix(i, j);
156 }
157
158 // Return the mass matrix.
ComputeMmatrixGlobal(ChMatrixRef M)159 void ChElementShellANCF_3423::ComputeMmatrixGlobal(ChMatrixRef M) {
160 M = m_MassMatrix;
161 }
162
163 // -----------------------------------------------------------------------------
164 // Mass matrix calculation
165 // -----------------------------------------------------------------------------
166
167 /// This class defines the calculations for the integrand of the inertia matrix.
168 class ShellANCF_Mass : public ChIntegrable3D<ChMatrixNM<double, 24, 24>> {
169 public:
ShellANCF_Mass(ChElementShellANCF_3423 * element)170 ShellANCF_Mass(ChElementShellANCF_3423* element) : m_element(element) {}
~ShellANCF_Mass()171 ~ShellANCF_Mass() {}
172
173 private:
174 ChElementShellANCF_3423* m_element;
175
176 virtual void Evaluate(ChMatrixNM<double, 24, 24>& result, const double x, const double y, const double z) override;
177 };
178
Evaluate(ChMatrixNM<double,24,24> & result,const double x,const double y,const double z)179 void ShellANCF_Mass::Evaluate(ChMatrixNM<double, 24, 24>& result, const double x, const double y, const double z) {
180 ChElementShellANCF_3423::ShapeVector N;
181 m_element->ShapeFunctions(N, x, y, z);
182
183 // S=[N1*eye(3) N2*eye(3) N3*eye(3) N4*eye(3) N5*eye(3) N6*eye(3) N7*eye(3) N8*eye(3)]
184 ChMatrixNM<double, 3, 24> S;
185 ChMatrix33<> Si;
186 S.setZero();
187 for (int i = 0; i < 8; i++) {
188 S(0, 3 * i + 0) = N(i);
189 S(1, 3 * i + 1) = N(i);
190 S(2, 3 * i + 2) = N(i);
191 }
192
193 double detJ0 = m_element->Calc_detJ0(x, y, z);
194
195 // perform r = S'*S, scaled by integration weights
196 result = detJ0 * m_element->m_GaussScaling * S.transpose() * S;
197 };
198
ComputeMassMatrix()199 void ChElementShellANCF_3423::ComputeMassMatrix() {
200 m_MassMatrix.setZero();
201
202 for (size_t kl = 0; kl < m_numLayers; kl++) {
203 double rho = m_layers[kl].GetMaterial()->Get_rho();
204 ShellANCF_Mass myformula(this);
205 ChMatrixNM<double, 24, 24> TempMassMatrix;
206 TempMassMatrix.setZero();
207 ChQuadrature::Integrate3D<ChMatrixNM<double, 24, 24>>(TempMassMatrix, // result of integration will go there
208 myformula, // formula to integrate
209 -1, 1, // x limits
210 -1, 1, // y limits
211 m_GaussZ[kl], m_GaussZ[kl + 1], // z limits
212 2 // order of integration
213 );
214 TempMassMatrix *= rho;
215 m_MassMatrix += TempMassMatrix;
216 }
217 }
218 /// This class computes and adds corresponding masses to ElementGeneric member m_TotalMass
ComputeNodalMass()219 void ChElementShellANCF_3423::ComputeNodalMass() {
220 m_nodes[0]->m_TotalMass += m_MassMatrix(0, 0) + m_MassMatrix(0, 6) + m_MassMatrix(0, 12) + m_MassMatrix(0, 18);
221 m_nodes[1]->m_TotalMass += m_MassMatrix(6, 6) + m_MassMatrix(6, 0) + m_MassMatrix(6, 12) + m_MassMatrix(6, 18);
222 m_nodes[2]->m_TotalMass += m_MassMatrix(12, 12) + m_MassMatrix(12, 0) + m_MassMatrix(12, 6) + m_MassMatrix(12, 18);
223 m_nodes[3]->m_TotalMass += m_MassMatrix(18, 18) + m_MassMatrix(18, 0) + m_MassMatrix(18, 6) + m_MassMatrix(18, 12);
224 }
225 // -----------------------------------------------------------------------------
226 // Gravitational force calculation
227 // -----------------------------------------------------------------------------
228
229 /// This class defines the calculations for the integrand of the element gravity forces
230 class ShellANCF_Gravity : public ChIntegrable3D<ChElementShellANCF_3423::VectorN> {
231 public:
ShellANCF_Gravity(ChElementShellANCF_3423 * element)232 ShellANCF_Gravity(ChElementShellANCF_3423* element) : m_element(element) {}
~ShellANCF_Gravity()233 ~ShellANCF_Gravity() {}
234
235 private:
236 ChElementShellANCF_3423* m_element;
237
238 virtual void Evaluate(ChElementShellANCF_3423::VectorN& result,
239 const double x,
240 const double y,
241 const double z) override;
242 };
243
Evaluate(ChElementShellANCF_3423::VectorN & result,const double x,const double y,const double z)244 void ShellANCF_Gravity::Evaluate(ChElementShellANCF_3423::VectorN& result,
245 const double x,
246 const double y,
247 const double z) {
248 ChElementShellANCF_3423::ShapeVector N;
249 m_element->ShapeFunctions(N, x, y, z);
250
251 result = m_element->Calc_detJ0(x, y, z) * m_element->m_GaussScaling * N.transpose();
252 };
253
ComputeGravityForceScale()254 void ChElementShellANCF_3423::ComputeGravityForceScale() {
255 m_GravForceScale.setZero();
256
257 for (size_t kl = 0; kl < m_numLayers; kl++) {
258 double rho = m_layers[kl].GetMaterial()->Get_rho();
259 ShellANCF_Gravity myformula(this);
260 VectorN Fgravity;
261 Fgravity.setZero();
262 ChQuadrature::Integrate3D<VectorN>(Fgravity, // result of integration will go there
263 myformula, // formula to integrate
264 -1, 1, // x limits
265 -1, 1, // y limits
266 m_GaussZ[kl], m_GaussZ[kl + 1], // z limits
267 2 // order of integration
268 );
269
270 m_GravForceScale += rho * Fgravity;
271 }
272 }
273
274 // Compute the generalized force vector due to gravity using the efficient ANCF specific method
ComputeGravityForces(ChVectorDynamic<> & Fg,const ChVector<> & G_acc)275 void ChElementShellANCF_3423::ComputeGravityForces(ChVectorDynamic<>& Fg, const ChVector<>& G_acc) {
276 assert(Fg.size() == 3 * NSF);
277
278 // Calculate and add the generalized force due to gravity to the generalized internal force vector for the element.
279 // The generalized force due to gravity could be computed once prior to the start of the simulation if gravity was
280 // assumed constant throughout the entire simulation. However, this implementation assumes that the acceleration
281 // due to gravity, while a constant for the entire system, can change from step to step which could be useful for
282 // gravity loaded units tests as an example. The generalized force due to gravity is calculated in compact matrix
283 // form and is pre-mapped to the desired vector format
284 Eigen::Map<MatrixNx3> GravForceCompact(Fg.data(), NSF, 3);
285 GravForceCompact = m_GravForceScale * G_acc.eigen().transpose();
286 }
287
288 // -----------------------------------------------------------------------------
289 // Elastic force calculation
290 // -----------------------------------------------------------------------------
291
292 // The class ShellANCF_Force provides the integrand for the calculation of the internal forces
293 // for one layer of an ANCF shell element.
294 // The first 24 entries in the integrand represent the internal force.
295 // The next 5 entries represent the residual of the EAS nonlinear system.
296 // The last 25 entries represent the 5x5 Jacobian of the EAS nonlinear system.
297 // Capabilities of this class include: application of enhanced assumed strain (EAS) and
298 // assumed natural strain (ANS) formulations to avoid thickness and (transverse and in-plane)
299 // shear locking. This implementation also features a composite material implementation
300 // that allows for selecting a number of layers over the element thickness; each of which
301 // has an independent, user-selected fiber angle (direction for orthotropic constitutive behavior)
302 class ShellANCF_Force : public ChIntegrable3D<ChVectorN<double, 54>> {
303 public:
ShellANCF_Force(ChElementShellANCF_3423 * element,size_t kl,ChVectorN<double,5> * alpha_eas)304 ShellANCF_Force(ChElementShellANCF_3423* element, // Containing element
305 size_t kl, // Current layer index
306 ChVectorN<double, 5>* alpha_eas // Vector of internal parameters for EAS formulation
307 )
308 : m_element(element), m_kl(kl), m_alpha_eas(alpha_eas) {}
~ShellANCF_Force()309 ~ShellANCF_Force() {}
310
311 private:
312 ChElementShellANCF_3423* m_element;
313 size_t m_kl;
314 ChVectorN<double, 5>* m_alpha_eas;
315
316 /// Evaluate (strainD'*strain) at point x, include ANS and EAS.
317 virtual void Evaluate(ChVectorN<double, 54>& result, const double x, const double y, const double z) override;
318 };
319
Evaluate(ChVectorN<double,54> & result,const double x,const double y,const double z)320 void ShellANCF_Force::Evaluate(ChVectorN<double, 54>& result, const double x, const double y, const double z) {
321 // Element shape function
322 ChElementShellANCF_3423::ShapeVector N;
323 m_element->ShapeFunctions(N, x, y, z);
324
325 // Determinant of position vector gradient matrix: Initial configuration
326 ChElementShellANCF_3423::ShapeVector Nx;
327 ChElementShellANCF_3423::ShapeVector Ny;
328 ChElementShellANCF_3423::ShapeVector Nz;
329 ChMatrixNM<double, 1, 3> Nx_d0;
330 ChMatrixNM<double, 1, 3> Ny_d0;
331 ChMatrixNM<double, 1, 3> Nz_d0;
332 double detJ0 = m_element->Calc_detJ0(x, y, z, Nx, Ny, Nz, Nx_d0, Ny_d0, Nz_d0);
333
334 // ANS shape function
335 ChMatrixNM<double, 1, 4> S_ANS; // Shape function vector for Assumed Natural Strain
336 ChMatrixNM<double, 6, 5> M; // Shape function vector for Enhanced Assumed Strain
337 m_element->ShapeFunctionANSbilinearShell(S_ANS, x, y);
338 m_element->Basis_M(M, x, y, z);
339
340 // Transformation : Orthogonal transformation (A and J)
341 ChVector<double> G1xG2; // Cross product of first and second column of
342 double G1dotG1; // Dot product of first column of position vector gradient
343
344 G1xG2.x() = Nx_d0(1) * Ny_d0(2) - Nx_d0(2) * Ny_d0(1);
345 G1xG2.y() = Nx_d0(2) * Ny_d0(0) - Nx_d0(0) * Ny_d0(2);
346 G1xG2.z() = Nx_d0(0) * Ny_d0(1) - Nx_d0(1) * Ny_d0(0);
347 G1dotG1 = Nx_d0(0) * Nx_d0(0) + Nx_d0(1) * Nx_d0(1) + Nx_d0(2) * Nx_d0(2);
348
349 // Tangent Frame
350 ChVector<double> A1;
351 ChVector<double> A2;
352 ChVector<double> A3;
353 A1.x() = Nx_d0(0);
354 A1.y() = Nx_d0(1);
355 A1.z() = Nx_d0(2);
356 A1 = A1 / sqrt(G1dotG1);
357 A3 = G1xG2.GetNormalized();
358 A2.Cross(A3, A1);
359
360 // Direction for orthotropic material
361 double theta = m_element->GetLayer(m_kl).Get_theta(); // Fiber angle
362 ChVector<double> AA1;
363 ChVector<double> AA2;
364 ChVector<double> AA3;
365 AA1 = A1 * cos(theta) + A2 * sin(theta);
366 AA2 = -A1 * sin(theta) + A2 * cos(theta);
367 AA3 = A3;
368
369 /// Beta
370 ChMatrixNM<double, 3, 3> j0;
371 ChVector<double> j01;
372 ChVector<double> j02;
373 ChVector<double> j03;
374 ChVectorN<double, 9> beta;
375 // Calculates inverse of rd0 (j0) (position vector gradient: Initial Configuration)
376 j0(0, 0) = Ny_d0(1) * Nz_d0(2) - Nz_d0(1) * Ny_d0(2);
377 j0(0, 1) = Ny_d0(2) * Nz_d0(0) - Ny_d0(0) * Nz_d0(2);
378 j0(0, 2) = Ny_d0(0) * Nz_d0(1) - Nz_d0(0) * Ny_d0(1);
379 j0(1, 0) = Nz_d0(1) * Nx_d0(2) - Nx_d0(1) * Nz_d0(2);
380 j0(1, 1) = Nz_d0(2) * Nx_d0(0) - Nx_d0(2) * Nz_d0(0);
381 j0(1, 2) = Nz_d0(0) * Nx_d0(1) - Nz_d0(1) * Nx_d0(0);
382 j0(2, 0) = Nx_d0(1) * Ny_d0(2) - Ny_d0(1) * Nx_d0(2);
383 j0(2, 1) = Ny_d0(0) * Nx_d0(2) - Nx_d0(0) * Ny_d0(2);
384 j0(2, 2) = Nx_d0(0) * Ny_d0(1) - Ny_d0(0) * Nx_d0(1);
385 j0 /= detJ0;
386
387 j01[0] = j0(0, 0);
388 j02[0] = j0(1, 0);
389 j03[0] = j0(2, 0);
390 j01[1] = j0(0, 1);
391 j02[1] = j0(1, 1);
392 j03[1] = j0(2, 1);
393 j01[2] = j0(0, 2);
394 j02[2] = j0(1, 2);
395 j03[2] = j0(2, 2);
396
397 // Coefficients of contravariant transformation
398 beta(0) = Vdot(AA1, j01);
399 beta(1) = Vdot(AA2, j01);
400 beta(2) = Vdot(AA3, j01);
401 beta(3) = Vdot(AA1, j02);
402 beta(4) = Vdot(AA2, j02);
403 beta(5) = Vdot(AA3, j02);
404 beta(6) = Vdot(AA1, j03);
405 beta(7) = Vdot(AA2, j03);
406 beta(8) = Vdot(AA3, j03);
407
408 // Transformation matrix, function of fiber angle
409 const ChMatrixNM<double, 6, 6>& T0 = m_element->GetLayer(m_kl).Get_T0();
410 // Determinant of the initial position vector gradient at the element center
411 double detJ0C = m_element->GetLayer(m_kl).Get_detJ0C();
412
413 // Enhanced Assumed Strain
414 ChMatrixNM<double, 6, 5> G = T0 * M * (detJ0C / detJ0);
415 ChVectorN<double, 6> strain_EAS = G * (*m_alpha_eas);
416
417 ChVectorN<double, 8> ddNx = m_element->m_ddT * Nx.transpose();
418 ChVectorN<double, 8> ddNy = m_element->m_ddT * Ny.transpose();
419
420 ChVectorN<double, 8> d0d0Nx = m_element->m_d0d0T * Nx.transpose();
421 ChVectorN<double, 8> d0d0Ny = m_element->m_d0d0T * Ny.transpose();
422
423 // Strain component
424 ChVectorN<double, 6> strain_til;
425 strain_til(0) = 0.5 * ((Nx * ddNx)(0, 0) - (Nx * d0d0Nx)(0, 0));
426 strain_til(1) = 0.5 * ((Ny * ddNy)(0, 0) - (Ny * d0d0Ny)(0, 0));
427 strain_til(2) = (Nx * ddNy)(0, 0) - (Nx * d0d0Ny)(0, 0);
428 strain_til(3) = N(0) * m_element->m_strainANS(0) + N(2) * m_element->m_strainANS(1) +
429 N(4) * m_element->m_strainANS(2) + N(6) * m_element->m_strainANS(3);
430 strain_til(4) = S_ANS(0, 2) * m_element->m_strainANS(6) + S_ANS(0, 3) * m_element->m_strainANS(7);
431 strain_til(5) = S_ANS(0, 0) * m_element->m_strainANS(4) + S_ANS(0, 1) * m_element->m_strainANS(5);
432
433 // For orthotropic material
434 ChVectorN<double, 6> strain;
435
436 strain(0) = strain_til(0) * beta(0) * beta(0) + strain_til(1) * beta(3) * beta(3) +
437 strain_til(2) * beta(0) * beta(3) + strain_til(3) * beta(6) * beta(6) +
438 strain_til(4) * beta(0) * beta(6) + strain_til(5) * beta(3) * beta(6);
439 strain(1) = strain_til(0) * beta(1) * beta(1) + strain_til(1) * beta(4) * beta(4) +
440 strain_til(2) * beta(1) * beta(4) + strain_til(3) * beta(7) * beta(7) +
441 strain_til(4) * beta(1) * beta(7) + strain_til(5) * beta(4) * beta(7);
442 strain(2) = strain_til(0) * 2.0 * beta(0) * beta(1) + strain_til(1) * 2.0 * beta(3) * beta(4) +
443 strain_til(2) * (beta(1) * beta(3) + beta(0) * beta(4)) + strain_til(3) * 2.0 * beta(6) * beta(7) +
444 strain_til(4) * (beta(1) * beta(6) + beta(0) * beta(7)) +
445 strain_til(5) * (beta(4) * beta(6) + beta(3) * beta(7));
446 strain(3) = strain_til(0) * beta(2) * beta(2) + strain_til(1) * beta(5) * beta(5) +
447 strain_til(2) * beta(2) * beta(5) + strain_til(3) * beta(8) * beta(8) +
448 strain_til(4) * beta(2) * beta(8) + strain_til(5) * beta(5) * beta(8);
449 strain(4) = strain_til(0) * 2.0 * beta(0) * beta(2) + strain_til(1) * 2.0 * beta(3) * beta(5) +
450 strain_til(2) * (beta(2) * beta(3) + beta(0) * beta(5)) + strain_til(3) * 2.0 * beta(6) * beta(8) +
451 strain_til(4) * (beta(2) * beta(6) + beta(0) * beta(8)) +
452 strain_til(5) * (beta(5) * beta(6) + beta(3) * beta(8));
453 strain(5) = strain_til(0) * 2.0 * beta(1) * beta(2) + strain_til(1) * 2.0 * beta(4) * beta(5) +
454 strain_til(2) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3) * 2.0 * beta(7) * beta(8) +
455 strain_til(4) * (beta(2) * beta(7) + beta(1) * beta(8)) +
456 strain_til(5) * (beta(5) * beta(7) + beta(4) * beta(8));
457
458 // Strain derivative component
459
460 ChMatrixNM<double, 6, 24> strainD_til;
461
462 ChMatrixNM<double, 1, 24> tempB;
463 ChMatrixNM<double, 1, 3> tempB3;
464 ChMatrixNM<double, 1, 3> tempB31;
465
466 tempB3 = Nx * m_element->m_d;
467 for (int i = 0; i < 8; i++) {
468 for (int j = 0; j < 3; j++) {
469 tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i);
470 }
471 }
472 strainD_til.row(0) = tempB;
473
474 tempB3 = Ny * m_element->m_d;
475 for (int i = 0; i < 8; i++) {
476 for (int j = 0; j < 3; j++) {
477 tempB(0, i * 3 + j) = tempB3(0, j) * Ny(0, i);
478 }
479 }
480 strainD_til.row(1) = tempB;
481
482 tempB31 = Nx * m_element->m_d;
483 for (int i = 0; i < 8; i++) {
484 for (int j = 0; j < 3; j++) {
485 tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i) + tempB31(0, j) * Ny(0, i);
486 }
487 }
488 strainD_til.row(2) = tempB;
489
490 tempB.setZero();
491 for (int i = 0; i < 4; i++) {
492 tempB += N(i * 2) * m_element->m_strainANS_D.row(i);
493 }
494 strainD_til.row(3) = tempB; // strainD for zz
495
496 tempB.setZero();
497 for (int i = 0; i < 2; i++) {
498 tempB += S_ANS(0, i + 2) * m_element->m_strainANS_D.row(i + 6);
499 }
500 strainD_til.row(4) = tempB; // strainD for xz
501
502 tempB.setZero();
503 for (int i = 0; i < 2; i++) {
504 tempB += S_ANS(0, i) * m_element->m_strainANS_D.row(i + 4);
505 }
506 strainD_til.row(5) = tempB; // strainD for yz
507
508 // For orthotropic material
509 ChMatrixNM<double, 6, 24> strainD; // Derivative of the strains w.r.t. the coordinates. Includes orthotropy
510 for (int ii = 0; ii < 24; ii++) {
511 strainD(0, ii) = strainD_til(0, ii) * beta(0) * beta(0) + strainD_til(1, ii) * beta(3) * beta(3) +
512 strainD_til(2, ii) * beta(0) * beta(3) + strainD_til(3, ii) * beta(6) * beta(6) +
513 strainD_til(4, ii) * beta(0) * beta(6) + strainD_til(5, ii) * beta(3) * beta(6);
514 strainD(1, ii) = strainD_til(0, ii) * beta(1) * beta(1) + strainD_til(1, ii) * beta(4) * beta(4) +
515 strainD_til(2, ii) * beta(1) * beta(4) + strainD_til(3, ii) * beta(7) * beta(7) +
516 strainD_til(4, ii) * beta(1) * beta(7) + strainD_til(5, ii) * beta(4) * beta(7);
517 strainD(2, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(1) + strainD_til(1, ii) * 2.0 * beta(3) * beta(4) +
518 strainD_til(2, ii) * (beta(1) * beta(3) + beta(0) * beta(4)) +
519 strainD_til(3, ii) * 2.0 * beta(6) * beta(7) +
520 strainD_til(4, ii) * (beta(1) * beta(6) + beta(0) * beta(7)) +
521 strainD_til(5, ii) * (beta(4) * beta(6) + beta(3) * beta(7));
522 strainD(3, ii) = strainD_til(0, ii) * beta(2) * beta(2) + strainD_til(1, ii) * beta(5) * beta(5) +
523 strainD_til(2, ii) * beta(2) * beta(5) + strainD_til(3, ii) * beta(8) * beta(8) +
524 strainD_til(4, ii) * beta(2) * beta(8) + strainD_til(5) * beta(5) * beta(8);
525 strainD(4, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(2) + strainD_til(1, ii) * 2.0 * beta(3) * beta(5) +
526 strainD_til(2, ii) * (beta(2) * beta(3) + beta(0) * beta(5)) +
527 strainD_til(3, ii) * 2.0 * beta(6) * beta(8) +
528 strainD_til(4, ii) * (beta(2) * beta(6) + beta(0) * beta(8)) +
529 strainD_til(5, ii) * (beta(5) * beta(6) + beta(3) * beta(8));
530 strainD(5, ii) = strainD_til(0, ii) * 2.0 * beta(1) * beta(2) + strainD_til(1, ii) * 2.0 * beta(4) * beta(5) +
531 strainD_til(2, ii) * (beta(2) * beta(4) + beta(1) * beta(5)) +
532 strainD_til(3, ii) * 2.0 * beta(7) * beta(8) +
533 strainD_til(4, ii) * (beta(2) * beta(7) + beta(1) * beta(8)) +
534 strainD_til(5, ii) * (beta(5) * beta(7) + beta(4) * beta(8));
535 }
536
537 // Enhanced Assumed Strain 2nd
538 strain += strain_EAS;
539
540 // Strain time derivative for structural damping
541 ChVectorN<double, 6> DEPS;
542 DEPS.setZero();
543 for (int ii = 0; ii < 24; ii++) {
544 DEPS(0) += strainD(0, ii) * m_element->m_d_dt(ii);
545 DEPS(1) += strainD(1, ii) * m_element->m_d_dt(ii);
546 DEPS(2) += strainD(2, ii) * m_element->m_d_dt(ii);
547 DEPS(3) += strainD(3, ii) * m_element->m_d_dt(ii);
548 DEPS(4) += strainD(4, ii) * m_element->m_d_dt(ii);
549 DEPS(5) += strainD(5, ii) * m_element->m_d_dt(ii);
550 }
551
552 // Add structural damping
553 strain += DEPS * m_element->m_Alpha;
554
555 // Matrix of elastic coefficients: the input assumes the material *could* be orthotropic
556 const ChMatrixNM<double, 6, 6>& E_eps = m_element->GetLayer(m_kl).GetMaterial()->Get_E_eps();
557
558 // Internal force calculation
559 ChVectorN<double, 24> Fint = (strainD.transpose() * E_eps * strain) * (detJ0 * m_element->m_GaussScaling);
560
561 // EAS terms
562 ChMatrixNM<double, 5, 6> temp56 = G.transpose() * E_eps;
563 ChVectorN<double, 5> HE = (temp56 * strain) * (detJ0 * m_element->m_GaussScaling); // EAS residual
564 ChMatrixNM<double, 5, 5> KALPHA = (temp56 * G) * (detJ0 * m_element->m_GaussScaling); // EAS Jacobian
565
566 /// Total result vector
567 result.segment(0, 24) = Fint;
568 result.segment(24, 5) = HE;
569 result.segment(29, 5 * 5) = Eigen::Map<ChVectorN<double, 5 * 5>>(KALPHA.data(), 5 * 5);
570 }
571
ComputeInternalForces(ChVectorDynamic<> & Fi)572 void ChElementShellANCF_3423::ComputeInternalForces(ChVectorDynamic<>& Fi) {
573 // Current nodal coordinates and velocities
574 CalcCoordMatrix(m_d);
575 CalcCoordDerivMatrix(m_d_dt);
576 m_ddT = m_d * m_d.transpose();
577 // Assumed Natural Strain (ANS): Calculate m_strainANS and m_strainANS_D
578 CalcStrainANSbilinearShell();
579
580 Fi.setZero();
581
582 for (size_t kl = 0; kl < m_numLayers; kl++) {
583 ChVectorN<double, 24> Finternal;
584 ChVectorN<double, 5> HE;
585 ChMatrixNM<double, 5, 5> KALPHA;
586
587 // Initial guess for EAS parameters
588 ChVectorN<double, 5> alphaEAS = m_alphaEAS[kl];
589
590 // Newton loop for EAS
591 for (int count = 0; count < m_maxIterationsEAS; count++) {
592 ShellANCF_Force formula(this, kl, &alphaEAS);
593 ChVectorN<double, 54> result;
594 result.setZero();
595 ChQuadrature::Integrate3D<ChVectorN<double, 54>>(result, // result of integration
596 formula, // integrand formula
597 -1, 1, // x limits
598 -1, 1, // y limits
599 m_GaussZ[kl], m_GaussZ[kl + 1], // z limits
600 2 // order of integration
601 );
602
603 // Extract vectors and matrices from result of integration
604 Finternal = result.segment(0, 24);
605 HE = result.segment(24, 5);
606 KALPHA = Eigen::Map<ChMatrixNM<double, 5, 5>>(result.segment(29, 5 * 5).data(), 5, 5);
607
608 // Check convergence (residual check)
609 double norm_HE = HE.norm();
610 if (norm_HE < m_toleranceEAS)
611 break;
612
613 // Calculate increment and update EAS parameters
614 ChVectorN<double, 5> sol = KALPHA.colPivHouseholderQr().solve(HE);
615 alphaEAS -= sol;
616
617 if (count >= 2)
618 GetLog() << " count " << count << " NormHE " << norm_HE << "\n";
619 }
620
621 // Accumulate internal force
622 Fi -= Finternal;
623
624 // Cache alphaEAS and KALPHA for use in Jacobian calculation
625 m_alphaEAS[kl] = alphaEAS;
626 m_KalphaEAS[kl] = KALPHA;
627
628 } // Layer Loop
629 }
630
631 // -----------------------------------------------------------------------------
632 // Jacobians of internal forces
633 // -----------------------------------------------------------------------------
634
635 // The class ShellANCF_Jacobian provides the integrand for the calculation of the Jacobians
636 // (stiffness and damping matrices) of the internal forces for one layer of an ANCF
637 // shell element.
638 // The first 576 entries in the integrated vector represent the 24x24 Jacobian
639 // Kfactor * [K] + Rfactor * [R]
640 // where K does not include the EAS contribution.
641 // The last 120 entries represent the 5x24 cross-dependency matrix.
642 class ShellANCF_Jacobian : public ChIntegrable3D<ChVectorN<double, 696>> {
643 public:
ShellANCF_Jacobian(ChElementShellANCF_3423 * element,double Kfactor,double Rfactor,size_t kl)644 ShellANCF_Jacobian(ChElementShellANCF_3423* element, // Containing element
645 double Kfactor, // Scaling coefficient for stiffness component
646 double Rfactor, // Scaling coefficient for damping component
647 size_t kl // Current layer index
648 )
649 : m_element(element), m_Kfactor(Kfactor), m_Rfactor(Rfactor), m_kl(kl) {}
650
651 private:
652 ChElementShellANCF_3423* m_element;
653 double m_Kfactor;
654 double m_Rfactor;
655 size_t m_kl;
656
657 // Evaluate integrand at the specified point.
658 virtual void Evaluate(ChVectorN<double, 696>& result, const double x, const double y, const double z) override;
659 };
660
Evaluate(ChVectorN<double,696> & result,const double x,const double y,const double z)661 void ShellANCF_Jacobian::Evaluate(ChVectorN<double, 696>& result, const double x, const double y, const double z) {
662 // Element shape function
663 ChElementShellANCF_3423::ShapeVector N;
664 m_element->ShapeFunctions(N, x, y, z);
665
666 // Determinant of position vector gradient matrix: Initial configuration
667 ChElementShellANCF_3423::ShapeVector Nx;
668 ChElementShellANCF_3423::ShapeVector Ny;
669 ChElementShellANCF_3423::ShapeVector Nz;
670 ChMatrixNM<double, 1, 3> Nx_d0;
671 ChMatrixNM<double, 1, 3> Ny_d0;
672 ChMatrixNM<double, 1, 3> Nz_d0;
673 double detJ0 = m_element->Calc_detJ0(x, y, z, Nx, Ny, Nz, Nx_d0, Ny_d0, Nz_d0);
674
675 // ANS shape function
676 ChMatrixNM<double, 1, 4> S_ANS; // Shape function vector for Assumed Natural Strain
677 ChMatrixNM<double, 6, 5> M; // Shape function vector for Enhanced Assumed Strain
678 m_element->ShapeFunctionANSbilinearShell(S_ANS, x, y);
679 m_element->Basis_M(M, x, y, z);
680
681 // Transformation : Orthogonal transformation (A and J)
682 ChVector<double> G1xG2; // Cross product of first and second column of
683 double G1dotG1; // Dot product of first column of position vector gradient
684
685 G1xG2.x() = Nx_d0(1) * Ny_d0(2) - Nx_d0(2) * Ny_d0(1);
686 G1xG2.y() = Nx_d0(2) * Ny_d0(0) - Nx_d0(0) * Ny_d0(2);
687 G1xG2.z() = Nx_d0(0) * Ny_d0(1) - Nx_d0(1) * Ny_d0(0);
688 G1dotG1 = Nx_d0(0) * Nx_d0(0) + Nx_d0(1) * Nx_d0(1) + Nx_d0(2) * Nx_d0(2);
689
690 // Tangent Frame
691 ChVector<double> A1;
692 ChVector<double> A2;
693 ChVector<double> A3;
694 A1.x() = Nx_d0(0);
695 A1.y() = Nx_d0(1);
696 A1.z() = Nx_d0(2);
697 A1 = A1 / sqrt(G1dotG1);
698 A3 = G1xG2.GetNormalized();
699 A2.Cross(A3, A1);
700
701 // Direction for orthotropic material
702 double theta = m_element->GetLayer(m_kl).Get_theta(); // Fiber angle
703 ChVector<double> AA1;
704 ChVector<double> AA2;
705 ChVector<double> AA3;
706 AA1 = A1 * cos(theta) + A2 * sin(theta);
707 AA2 = -A1 * sin(theta) + A2 * cos(theta);
708 AA3 = A3;
709
710 /// Beta
711 ChMatrixNM<double, 3, 3> j0;
712 ChVector<double> j01;
713 ChVector<double> j02;
714 ChVector<double> j03;
715 // Calculates inverse of rd0 (j0) (position vector gradient: Initial Configuration)
716 j0(0, 0) = Ny_d0(1) * Nz_d0(2) - Nz_d0(1) * Ny_d0(2);
717 j0(0, 1) = Ny_d0(2) * Nz_d0(0) - Ny_d0(0) * Nz_d0(2);
718 j0(0, 2) = Ny_d0(0) * Nz_d0(1) - Nz_d0(0) * Ny_d0(1);
719 j0(1, 0) = Nz_d0(1) * Nx_d0(2) - Nx_d0(1) * Nz_d0(2);
720 j0(1, 1) = Nz_d0(2) * Nx_d0(0) - Nx_d0(2) * Nz_d0(0);
721 j0(1, 2) = Nz_d0(0) * Nx_d0(1) - Nz_d0(1) * Nx_d0(0);
722 j0(2, 0) = Nx_d0(1) * Ny_d0(2) - Ny_d0(1) * Nx_d0(2);
723 j0(2, 1) = Ny_d0(0) * Nx_d0(2) - Nx_d0(0) * Ny_d0(2);
724 j0(2, 2) = Nx_d0(0) * Ny_d0(1) - Ny_d0(0) * Nx_d0(1);
725 j0 /= detJ0;
726
727 j01[0] = j0(0, 0);
728 j02[0] = j0(1, 0);
729 j03[0] = j0(2, 0);
730 j01[1] = j0(0, 1);
731 j02[1] = j0(1, 1);
732 j03[1] = j0(2, 1);
733 j01[2] = j0(0, 2);
734 j02[2] = j0(1, 2);
735 j03[2] = j0(2, 2);
736
737 // Coefficients of contravariant transformation
738 ChVectorN<double, 9> beta;
739 beta(0) = Vdot(AA1, j01);
740 beta(1) = Vdot(AA2, j01);
741 beta(2) = Vdot(AA3, j01);
742 beta(3) = Vdot(AA1, j02);
743 beta(4) = Vdot(AA2, j02);
744 beta(5) = Vdot(AA3, j02);
745 beta(6) = Vdot(AA1, j03);
746 beta(7) = Vdot(AA2, j03);
747 beta(8) = Vdot(AA3, j03);
748
749 // Transformation matrix, function of fiber angle
750 const ChMatrixNM<double, 6, 6>& T0 = m_element->GetLayer(m_kl).Get_T0();
751 // Determinant of the initial position vector gradient at the element center
752 double detJ0C = m_element->GetLayer(m_kl).Get_detJ0C();
753
754 // Enhanced Assumed Strain
755 ChMatrixNM<double, 6, 5> G = T0 * M * (detJ0C / detJ0);
756 ChVectorN<double, 6> strain_EAS = G * m_element->m_alphaEAS[m_kl];
757
758 ChVectorN<double, 8> ddNx = m_element->m_ddT * Nx.transpose();
759 ChVectorN<double, 8> ddNy = m_element->m_ddT * Ny.transpose();
760
761 ChVectorN<double, 8> d0d0Nx = m_element->m_d0d0T * Nx.transpose();
762 ChVectorN<double, 8> d0d0Ny = m_element->m_d0d0T * Ny.transpose();
763
764 // Strain component
765 ChVectorN<double, 6> strain_til;
766 strain_til(0) = 0.5 * ((Nx * ddNx)(0, 0) - (Nx * d0d0Nx)(0, 0));
767 strain_til(1) = 0.5 * ((Ny * ddNy)(0, 0) - (Ny * d0d0Ny)(0, 0));
768 strain_til(2) = (Nx * ddNy)(0, 0) - (Nx * d0d0Ny)(0, 0);
769 strain_til(3) = N(0) * m_element->m_strainANS(0) + N(2) * m_element->m_strainANS(1) +
770 N(4) * m_element->m_strainANS(2) + N(6) * m_element->m_strainANS(3);
771 strain_til(4) = S_ANS(0, 2) * m_element->m_strainANS(6) + S_ANS(0, 3) * m_element->m_strainANS(7);
772 strain_til(5) = S_ANS(0, 0) * m_element->m_strainANS(4) + S_ANS(0, 1) * m_element->m_strainANS(5);
773
774 // For orthotropic material
775 ChVectorN<double, 6> strain;
776
777 strain(0) = strain_til(0) * beta(0) * beta(0) + strain_til(1) * beta(3) * beta(3) +
778 strain_til(2) * beta(0) * beta(3) + strain_til(3) * beta(6) * beta(6) +
779 strain_til(4) * beta(0) * beta(6) + strain_til(5) * beta(3) * beta(6);
780 strain(1) = strain_til(0) * beta(1) * beta(1) + strain_til(1) * beta(4) * beta(4) +
781 strain_til(2) * beta(1) * beta(4) + strain_til(3) * beta(7) * beta(7) +
782 strain_til(4) * beta(1) * beta(7) + strain_til(5) * beta(4) * beta(7);
783 strain(2) = strain_til(0) * 2.0 * beta(0) * beta(1) + strain_til(1) * 2.0 * beta(3) * beta(4) +
784 strain_til(2) * (beta(1) * beta(3) + beta(0) * beta(4)) + strain_til(3) * 2.0 * beta(6) * beta(7) +
785 strain_til(4) * (beta(1) * beta(6) + beta(0) * beta(7)) +
786 strain_til(5) * (beta(4) * beta(6) + beta(3) * beta(7));
787 strain(3) = strain_til(0) * beta(2) * beta(2) + strain_til(1) * beta(5) * beta(5) +
788 strain_til(2) * beta(2) * beta(5) + strain_til(3) * beta(8) * beta(8) +
789 strain_til(4) * beta(2) * beta(8) + strain_til(5) * beta(5) * beta(8);
790 strain(4) = strain_til(0) * 2.0 * beta(0) * beta(2) + strain_til(1) * 2.0 * beta(3) * beta(5) +
791 strain_til(2) * (beta(2) * beta(3) + beta(0) * beta(5)) + strain_til(3) * 2.0 * beta(6) * beta(8) +
792 strain_til(4) * (beta(2) * beta(6) + beta(0) * beta(8)) +
793 strain_til(5) * (beta(5) * beta(6) + beta(3) * beta(8));
794 strain(5) = strain_til(0) * 2.0 * beta(1) * beta(2) + strain_til(1) * 2.0 * beta(4) * beta(5) +
795 strain_til(2) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3) * 2.0 * beta(7) * beta(8) +
796 strain_til(4) * (beta(2) * beta(7) + beta(1) * beta(8)) +
797 strain_til(5) * (beta(5) * beta(7) + beta(4) * beta(8));
798
799 // Strain derivative component
800
801 ChMatrixNM<double, 6, 24> strainD_til;
802 strainD_til.setZero();
803
804 ChMatrixNM<double, 1, 24> tempB;
805 ChMatrixNM<double, 1, 3> tempB3;
806 ChMatrixNM<double, 1, 3> tempB31;
807
808 tempB3 = Nx * m_element->m_d;
809 for (int i = 0; i < 8; i++) {
810 for (int j = 0; j < 3; j++) {
811 tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i);
812 }
813 }
814 strainD_til.row(0) = tempB;
815
816 tempB3 = Ny * m_element->m_d;
817 for (int i = 0; i < 8; i++) {
818 for (int j = 0; j < 3; j++) {
819 tempB(0, i * 3 + j) = tempB3(0, j) * Ny(0, i);
820 }
821 }
822 strainD_til.row(1) = tempB;
823
824 tempB31 = Nx * m_element->m_d;
825 for (int i = 0; i < 8; i++) {
826 for (int j = 0; j < 3; j++) {
827 tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i) + tempB31(0, j) * Ny(0, i);
828 }
829 }
830 strainD_til.row(2) = tempB;
831
832 tempB.setZero();
833 for (int i = 0; i < 4; i++) {
834 tempB += N(i * 2) * m_element->m_strainANS_D.row(i);
835 }
836 strainD_til.row(3) = tempB; // strainD for zz
837
838 tempB.setZero();
839 for (int i = 0; i < 2; i++) {
840 tempB += S_ANS(0, i + 2) * m_element->m_strainANS_D.row(i + 6);
841 }
842 strainD_til.row(4) = tempB; // strainD for xz
843
844 tempB.setZero();
845 for (int i = 0; i < 2; i++) {
846 tempB += S_ANS(0, i) * m_element->m_strainANS_D.row(i + 4);
847 }
848 strainD_til.row(5) = tempB; // strainD for yz
849
850 //// For orthotropic material
851 ChMatrixNM<double, 6, 24> strainD; // Derivative of the strains w.r.t. the coordinates. Includes orthotropy
852 for (int ii = 0; ii < 24; ii++) {
853 strainD(0, ii) = strainD_til(0, ii) * beta(0) * beta(0) + strainD_til(1, ii) * beta(3) * beta(3) +
854 strainD_til(2, ii) * beta(0) * beta(3) + strainD_til(3, ii) * beta(6) * beta(6) +
855 strainD_til(4, ii) * beta(0) * beta(6) + strainD_til(5, ii) * beta(3) * beta(6);
856 strainD(1, ii) = strainD_til(0, ii) * beta(1) * beta(1) + strainD_til(1, ii) * beta(4) * beta(4) +
857 strainD_til(2, ii) * beta(1) * beta(4) + strainD_til(3, ii) * beta(7) * beta(7) +
858 strainD_til(4, ii) * beta(1) * beta(7) + strainD_til(5, ii) * beta(4) * beta(7);
859 strainD(2, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(1) + strainD_til(1, ii) * 2.0 * beta(3) * beta(4) +
860 strainD_til(2, ii) * (beta(1) * beta(3) + beta(0) * beta(4)) +
861 strainD_til(3, ii) * 2.0 * beta(6) * beta(7) +
862 strainD_til(4, ii) * (beta(1) * beta(6) + beta(0) * beta(7)) +
863 strainD_til(5, ii) * (beta(4) * beta(6) + beta(3) * beta(7));
864 strainD(3, ii) = strainD_til(0, ii) * beta(2) * beta(2) + strainD_til(1, ii) * beta(5) * beta(5) +
865 strainD_til(2, ii) * beta(2) * beta(5) + strainD_til(3, ii) * beta(8) * beta(8) +
866 strainD_til(4, ii) * beta(2) * beta(8) + strainD_til(5) * beta(5) * beta(8);
867 strainD(4, ii) = strainD_til(0, ii) * 2.0 * beta(0) * beta(2) + strainD_til(1, ii) * 2.0 * beta(3) * beta(5) +
868 strainD_til(2, ii) * (beta(2) * beta(3) + beta(0) * beta(5)) +
869 strainD_til(3, ii) * 2.0 * beta(6) * beta(8) +
870 strainD_til(4, ii) * (beta(2) * beta(6) + beta(0) * beta(8)) +
871 strainD_til(5, ii) * (beta(5) * beta(6) + beta(3) * beta(8));
872 strainD(5, ii) = strainD_til(0, ii) * 2.0 * beta(1) * beta(2) + strainD_til(1, ii) * 2.0 * beta(4) * beta(5) +
873 strainD_til(2, ii) * (beta(2) * beta(4) + beta(1) * beta(5)) +
874 strainD_til(3, ii) * 2.0 * beta(7) * beta(8) +
875 strainD_til(4, ii) * (beta(2) * beta(7) + beta(1) * beta(8)) +
876 strainD_til(5, ii) * (beta(5) * beta(7) + beta(4) * beta(8));
877 }
878
879 /// Gd : Jacobian (w.r.t. coordinates) of the initial position vector gradient matrix
880 ChMatrixNM<double, 9, 24> Gd;
881 Gd.setZero();
882
883 for (int ii = 0; ii < 8; ii++) {
884 Gd(0, 3 * ii) = j0(0, 0) * Nx(0, ii) + j0(1, 0) * Ny(0, ii) + j0(2, 0) * Nz(0, ii);
885 Gd(1, 3 * ii + 1) = j0(0, 0) * Nx(0, ii) + j0(1, 0) * Ny(0, ii) + j0(2, 0) * Nz(0, ii);
886 Gd(2, 3 * ii + 2) = j0(0, 0) * Nx(0, ii) + j0(1, 0) * Ny(0, ii) + j0(2, 0) * Nz(0, ii);
887
888 Gd(3, 3 * ii) = j0(0, 1) * Nx(0, ii) + j0(1, 1) * Ny(0, ii) + j0(2, 1) * Nz(0, ii);
889 Gd(4, 3 * ii + 1) = j0(0, 1) * Nx(0, ii) + j0(1, 1) * Ny(0, ii) + j0(2, 1) * Nz(0, ii);
890 Gd(5, 3 * ii + 2) = j0(0, 1) * Nx(0, ii) + j0(1, 1) * Ny(0, ii) + j0(2, 1) * Nz(0, ii);
891
892 Gd(6, 3 * ii) = j0(0, 2) * Nx(0, ii) + j0(1, 2) * Ny(0, ii) + j0(2, 2) * Nz(0, ii);
893 Gd(7, 3 * ii + 1) = j0(0, 2) * Nx(0, ii) + j0(1, 2) * Ny(0, ii) + j0(2, 2) * Nz(0, ii);
894 Gd(8, 3 * ii + 2) = j0(0, 2) * Nx(0, ii) + j0(1, 2) * Ny(0, ii) + j0(2, 2) * Nz(0, ii);
895 }
896
897 // Enhanced Assumed Strain 2nd
898 strain += strain_EAS;
899
900 // Structural damping
901 // Strain time derivative for structural damping
902 ChVectorN<double, 6> DEPS;
903 DEPS.setZero();
904 for (int ii = 0; ii < 24; ii++) {
905 DEPS(0) += strainD(0, ii) * m_element->m_d_dt(ii);
906 DEPS(1) += strainD(1, ii) * m_element->m_d_dt(ii);
907 DEPS(2) += strainD(2, ii) * m_element->m_d_dt(ii);
908 DEPS(3) += strainD(3, ii) * m_element->m_d_dt(ii);
909 DEPS(4) += strainD(4, ii) * m_element->m_d_dt(ii);
910 DEPS(5) += strainD(5, ii) * m_element->m_d_dt(ii);
911 }
912
913 // Add structural damping
914 strain += DEPS * m_element->m_Alpha;
915
916 // Matrix of elastic coefficients: The input assumes the material *could* be orthotropic
917 const ChMatrixNM<double, 6, 6>& E_eps = m_element->GetLayer(m_kl).GetMaterial()->Get_E_eps();
918
919 // Stress tensor calculation
920 ChVectorN<double, 6> stress = E_eps * strain;
921
922 // Declaration and computation of Sigm, to be removed
923 ChMatrixNM<double, 9, 9> Sigm; ///< Rearrangement of stress vector (not always needed)
924 Sigm.setZero();
925
926 Sigm(0, 0) = stress(0); // XX
927 Sigm(1, 1) = stress(0);
928 Sigm(2, 2) = stress(0);
929
930 Sigm(0, 3) = stress(2); // XY
931 Sigm(1, 4) = stress(2);
932 Sigm(2, 5) = stress(2);
933
934 Sigm(0, 6) = stress(4); // XZ
935 Sigm(1, 7) = stress(4);
936 Sigm(2, 8) = stress(4);
937
938 Sigm(3, 0) = stress(2); // XY
939 Sigm(4, 1) = stress(2);
940 Sigm(5, 2) = stress(2);
941
942 Sigm(3, 3) = stress(1); // YY
943 Sigm(4, 4) = stress(1);
944 Sigm(5, 5) = stress(1);
945
946 Sigm(3, 6) = stress(5); // YZ
947 Sigm(4, 7) = stress(5);
948 Sigm(5, 8) = stress(5);
949
950 Sigm(6, 0) = stress(4); // XZ
951 Sigm(7, 1) = stress(4);
952 Sigm(8, 2) = stress(4);
953
954 Sigm(6, 3) = stress(5); // YZ
955 Sigm(7, 4) = stress(5);
956 Sigm(8, 5) = stress(5);
957
958 Sigm(6, 6) = stress(3); // ZZ
959 Sigm(7, 7) = stress(3);
960 Sigm(8, 8) = stress(3);
961
962 // Jacobian of internal forces (excluding the EAS contribution).
963 ChMatrixNM<double, 24, 24> KTE;
964 KTE = (strainD.transpose() * E_eps * strainD) * (m_Kfactor + m_Rfactor * m_element->m_Alpha) +
965 (Gd.transpose() * Sigm * Gd) * m_Kfactor;
966 KTE *= detJ0 * m_element->m_GaussScaling;
967
968 // EAS cross-dependency matrix.
969 ChMatrixNM<double, 5, 24> GDEPSP = (G.transpose() * E_eps * strainD) * (detJ0 * m_element->m_GaussScaling);
970
971 // Load result vector (integrand)
972 result.segment(0, 24 * 24) = Eigen::Map<ChVectorN<double, 24 * 24>>(KTE.data(), 24 * 24);
973 result.segment(576, 5 * 24) = Eigen::Map<ChVectorN<double, 5 * 24>>(GDEPSP.data(), 5 * 24);
974 }
975
ComputeInternalJacobians(double Kfactor,double Rfactor)976 void ChElementShellANCF_3423::ComputeInternalJacobians(double Kfactor, double Rfactor) {
977 // Note that the matrices with current nodal coordinates and velocities are
978 // already available in m_d and m_d_dt (as set in ComputeInternalForces).
979 // Similarly, the ANS strain and strain derivatives are already available in
980 // m_strainANS and m_strainANS_D (as calculated in ComputeInternalForces).
981
982 m_JacobianMatrix.setZero();
983
984 // Loop over all layers.
985 for (size_t kl = 0; kl < m_numLayers; kl++) {
986 ShellANCF_Jacobian formula(this, Kfactor, Rfactor, kl);
987 ChVectorN<double, 696> result;
988 result.setZero();
989 ChQuadrature::Integrate3D<ChVectorN<double, 696>>(result, // result of integration
990 formula, // integrand formula
991 -1, 1, // x limits
992 -1, 1, // y limits
993 m_GaussZ[kl], m_GaussZ[kl + 1], // z limits
994 2 // order of integration
995 );
996
997 // Extract matrices from result of integration
998 ChMatrixNM<double, 24, 24> KTE;
999 ChMatrixNM<double, 5, 24> GDEPSP;
1000 KTE = Eigen::Map<ChMatrixNM<double, 24, 24>>(result.segment(0, 24 * 24).data(), 24, 24);
1001 GDEPSP = Eigen::Map<ChMatrixNM<double, 5, 24>>(result.segment(576, 5 * 24).data(), 5, 24);
1002
1003 // Include EAS contribution to the stiffness component (hence scaled by Kfactor)
1004 // EAS = GDEPSP' * KalphaEAS_inv * GDEPSP
1005 ChMatrixNM<double, 5, 5> KalphaEAS_inv = m_KalphaEAS[kl].inverse();
1006 m_JacobianMatrix += KTE - Kfactor * GDEPSP.transpose() * KalphaEAS_inv * GDEPSP;
1007 }
1008 }
1009
1010 // -----------------------------------------------------------------------------
1011 // Shape functions
1012 // -----------------------------------------------------------------------------
1013
ShapeFunctions(ShapeVector & N,double x,double y,double z)1014 void ChElementShellANCF_3423::ShapeFunctions(ShapeVector& N, double x, double y, double z) {
1015 double c = m_thickness;
1016 N(0) = 0.25 * (1.0 - x) * (1.0 - y);
1017 N(1) = z * c / 2.0 * 0.25 * (1.0 - x) * (1.0 - y);
1018 N(2) = 0.25 * (1.0 + x) * (1.0 - y);
1019 N(3) = z * c / 2.0 * 0.25 * (1.0 + x) * (1.0 - y);
1020 N(4) = 0.25 * (1.0 + x) * (1.0 + y);
1021 N(5) = z * c / 2.0 * 0.25 * (1.0 + x) * (1.0 + y);
1022 N(6) = 0.25 * (1.0 - x) * (1.0 + y);
1023 N(7) = z * c / 2.0 * 0.25 * (1.0 - x) * (1.0 + y);
1024 }
1025
ShapeFunctionsDerivativeX(ShapeVector & Nx,double x,double y,double z)1026 void ChElementShellANCF_3423::ShapeFunctionsDerivativeX(ShapeVector& Nx, double x, double y, double z) {
1027 double a = GetLengthX();
1028 double c = m_thickness;
1029
1030 Nx(0) = 0.25 * (-2.0 / a) * (1.0 - y);
1031 Nx(1) = z * c / 2.0 * 0.25 * (-2.0 / a) * (1.0 - y);
1032 Nx(2) = 0.25 * (2.0 / a) * (1.0 - y);
1033 Nx(3) = z * c / 2.0 * 0.25 * (2.0 / a) * (1.0 - y);
1034 Nx(4) = 0.25 * (2.0 / a) * (1.0 + y);
1035 Nx(5) = z * c / 2.0 * 0.25 * (2.0 / a) * (1.0 + y);
1036 Nx(6) = 0.25 * (-2.0 / a) * (1.0 + y);
1037 Nx(7) = z * c / 2.0 * 0.25 * (-2.0 / a) * (1.0 + y);
1038 }
1039
ShapeFunctionsDerivativeY(ShapeVector & Ny,double x,double y,double z)1040 void ChElementShellANCF_3423::ShapeFunctionsDerivativeY(ShapeVector& Ny, double x, double y, double z) {
1041 double b = GetLengthY();
1042 double c = m_thickness;
1043
1044 Ny(0) = 0.25 * (1.0 - x) * (-2.0 / b);
1045 Ny(1) = z * c / 2.0 * 0.25 * (1.0 - x) * (-2.0 / b);
1046 Ny(2) = 0.25 * (1.0 + x) * (-2.0 / b);
1047 Ny(3) = z * c / 2.0 * 0.25 * (1.0 + x) * (-2.0 / b);
1048 Ny(4) = 0.25 * (1.0 + x) * (2.0 / b);
1049 Ny(5) = z * c / 2.0 * 0.25 * (1.0 + x) * (2.0 / b);
1050 Ny(6) = 0.25 * (1.0 - x) * (2.0 / b);
1051 Ny(7) = z * c / 2.0 * 0.25 * (1.0 - x) * (2.0 / b);
1052 }
1053
ShapeFunctionsDerivativeZ(ShapeVector & Nz,double x,double y,double z)1054 void ChElementShellANCF_3423::ShapeFunctionsDerivativeZ(ShapeVector& Nz, double x, double y, double z) {
1055 Nz(0) = 0.0;
1056 Nz(1) = 0.250 * (1.0 - x) * (1.0 - y);
1057 Nz(2) = 0.0;
1058 Nz(3) = 0.250 * (1.0 + x) * (1.0 - y);
1059 Nz(4) = 0.0;
1060 Nz(5) = 0.250 * (1.0 + x) * (1.0 + y);
1061 Nz(6) = 0.0;
1062 Nz(7) = 0.250 * (1.0 - x) * (1.0 + y);
1063 }
1064
Basis_M(ChMatrixNM<double,6,5> & M,double x,double y,double z)1065 void ChElementShellANCF_3423::Basis_M(ChMatrixNM<double, 6, 5>& M, double x, double y, double z) {
1066 M.setZero();
1067 M(0, 0) = x;
1068 M(1, 1) = y;
1069 M(2, 2) = x;
1070 M(2, 3) = y;
1071 M(3, 4) = z;
1072 }
1073
1074 // -----------------------------------------------------------------------------
1075 // Helper functions
1076 // -----------------------------------------------------------------------------
Calc_detJ0(double x,double y,double z,ShapeVector & Nx,ShapeVector & Ny,ShapeVector & Nz,ChMatrixNM<double,1,3> & Nx_d0,ChMatrixNM<double,1,3> & Ny_d0,ChMatrixNM<double,1,3> & Nz_d0)1077 double ChElementShellANCF_3423::Calc_detJ0(double x,
1078 double y,
1079 double z,
1080 ShapeVector& Nx,
1081 ShapeVector& Ny,
1082 ShapeVector& Nz,
1083 ChMatrixNM<double, 1, 3>& Nx_d0,
1084 ChMatrixNM<double, 1, 3>& Ny_d0,
1085 ChMatrixNM<double, 1, 3>& Nz_d0) {
1086 ShapeFunctionsDerivativeX(Nx, x, y, z);
1087 ShapeFunctionsDerivativeY(Ny, x, y, z);
1088 ShapeFunctionsDerivativeZ(Nz, x, y, z);
1089
1090 Nx_d0 = Nx * m_d0;
1091 Ny_d0 = Ny * m_d0;
1092 Nz_d0 = Nz * m_d0;
1093
1094 double detJ0 = Nx_d0(0, 0) * Ny_d0(0, 1) * Nz_d0(0, 2) + Ny_d0(0, 0) * Nz_d0(0, 1) * Nx_d0(0, 2) +
1095 Nz_d0(0, 0) * Nx_d0(0, 1) * Ny_d0(0, 2) - Nx_d0(0, 2) * Ny_d0(0, 1) * Nz_d0(0, 0) -
1096 Ny_d0(0, 2) * Nz_d0(0, 1) * Nx_d0(0, 0) - Nz_d0(0, 2) * Nx_d0(0, 1) * Ny_d0(0, 0);
1097
1098 return detJ0;
1099 }
1100
Calc_detJ0(double x,double y,double z)1101 double ChElementShellANCF_3423::Calc_detJ0(double x, double y, double z) {
1102 ChMatrixNM<double, 1, 8> Nx;
1103 ChMatrixNM<double, 1, 8> Ny;
1104 ChMatrixNM<double, 1, 8> Nz;
1105 ChMatrixNM<double, 1, 3> Nx_d0;
1106 ChMatrixNM<double, 1, 3> Ny_d0;
1107 ChMatrixNM<double, 1, 3> Nz_d0;
1108
1109 return Calc_detJ0(x, y, z, Nx, Ny, Nz, Nx_d0, Ny_d0, Nz_d0);
1110 }
1111
CalcCoordMatrix(ChMatrixNM<double,8,3> & d)1112 void ChElementShellANCF_3423::CalcCoordMatrix(ChMatrixNM<double, 8, 3>& d) {
1113 const ChVector<>& pA = m_nodes[0]->GetPos();
1114 const ChVector<>& dA = m_nodes[0]->GetD();
1115 const ChVector<>& pB = m_nodes[1]->GetPos();
1116 const ChVector<>& dB = m_nodes[1]->GetD();
1117 const ChVector<>& pC = m_nodes[2]->GetPos();
1118 const ChVector<>& dC = m_nodes[2]->GetD();
1119 const ChVector<>& pD = m_nodes[3]->GetPos();
1120 const ChVector<>& dD = m_nodes[3]->GetD();
1121
1122 d(0, 0) = pA.x();
1123 d(0, 1) = pA.y();
1124 d(0, 2) = pA.z();
1125 d(1, 0) = dA.x();
1126 d(1, 1) = dA.y();
1127 d(1, 2) = dA.z();
1128
1129 d(2, 0) = pB.x();
1130 d(2, 1) = pB.y();
1131 d(2, 2) = pB.z();
1132 d(3, 0) = dB.x();
1133 d(3, 1) = dB.y();
1134 d(3, 2) = dB.z();
1135
1136 d(4, 0) = pC.x();
1137 d(4, 1) = pC.y();
1138 d(4, 2) = pC.z();
1139 d(5, 0) = dC.x();
1140 d(5, 1) = dC.y();
1141 d(5, 2) = dC.z();
1142
1143 d(6, 0) = pD.x();
1144 d(6, 1) = pD.y();
1145 d(6, 2) = pD.z();
1146 d(7, 0) = dD.x();
1147 d(7, 1) = dD.y();
1148 d(7, 2) = dD.z();
1149 }
1150
CalcCoordDerivMatrix(ChVectorN<double,24> & dt)1151 void ChElementShellANCF_3423::CalcCoordDerivMatrix(ChVectorN<double, 24>& dt) {
1152 const ChVector<>& pA_dt = m_nodes[0]->GetPos_dt();
1153 const ChVector<>& dA_dt = m_nodes[0]->GetD_dt();
1154 const ChVector<>& pB_dt = m_nodes[1]->GetPos_dt();
1155 const ChVector<>& dB_dt = m_nodes[1]->GetD_dt();
1156 const ChVector<>& pC_dt = m_nodes[2]->GetPos_dt();
1157 const ChVector<>& dC_dt = m_nodes[2]->GetD_dt();
1158 const ChVector<>& pD_dt = m_nodes[3]->GetPos_dt();
1159 const ChVector<>& dD_dt = m_nodes[3]->GetD_dt();
1160
1161 dt(0) = pA_dt.x();
1162 dt(1) = pA_dt.y();
1163 dt(2) = pA_dt.z();
1164 dt(3) = dA_dt.x();
1165 dt(4) = dA_dt.y();
1166 dt(5) = dA_dt.z();
1167
1168 dt(6) = pB_dt.x();
1169 dt(7) = pB_dt.y();
1170 dt(8) = pB_dt.z();
1171 dt(9) = dB_dt.x();
1172 dt(10) = dB_dt.y();
1173 dt(11) = dB_dt.z();
1174
1175 dt(12) = pC_dt.x();
1176 dt(13) = pC_dt.y();
1177 dt(14) = pC_dt.z();
1178 dt(15) = dC_dt.x();
1179 dt(16) = dC_dt.y();
1180 dt(17) = dC_dt.z();
1181
1182 dt(18) = pD_dt.x();
1183 dt(19) = pD_dt.y();
1184 dt(20) = pD_dt.z();
1185 dt(21) = dD_dt.x();
1186 dt(22) = dD_dt.y();
1187 dt(23) = dD_dt.z();
1188 }
1189
1190 // -----------------------------------------------------------------------------
1191 // Assumed Natural Strain
1192 // -----------------------------------------------------------------------------
1193
1194 // ANS shape function (interpolation of strain and strainD in thickness direction)
ShapeFunctionANSbilinearShell(ChMatrixNM<double,1,4> & S_ANS,double x,double y)1195 void ChElementShellANCF_3423::ShapeFunctionANSbilinearShell(ChMatrixNM<double, 1, 4>& S_ANS, double x, double y) {
1196 S_ANS(0, 0) = -0.5 * x + 0.5;
1197 S_ANS(0, 1) = 0.5 * x + 0.5;
1198 S_ANS(0, 2) = -0.5 * y + 0.5;
1199 S_ANS(0, 3) = 0.5 * y + 0.5;
1200 }
1201
1202 // Calculate ANS strain and its Jacobian
CalcStrainANSbilinearShell()1203 void ChElementShellANCF_3423::CalcStrainANSbilinearShell() {
1204 std::vector<ChVector<>> knots(8);
1205
1206 knots[0] = ChVector<>(-1, -1, 0);
1207 knots[1] = ChVector<>(1, -1, 0);
1208 knots[2] = ChVector<>(-1, 1, 0);
1209 knots[3] = ChVector<>(1, 1, 0);
1210 knots[4] = ChVector<>(-1, 0, 0); // A
1211 knots[5] = ChVector<>(1, 0, 0); // B
1212 knots[6] = ChVector<>(0, -1, 0); // C
1213 knots[7] = ChVector<>(0, 1, 0); // D
1214
1215 ChMatrixNM<double, 1, 8> Nx;
1216 ChMatrixNM<double, 1, 8> Ny;
1217 ChMatrixNM<double, 1, 8> Nz;
1218 ChVectorN<double, 8> ddNz;
1219 ChVectorN<double, 8> d0d0Nz;
1220
1221 for (int kk = 0; kk < 8; kk++) {
1222 ShapeFunctionsDerivativeX(Nx, knots[kk].x(), knots[kk].y(), knots[kk].z());
1223 ShapeFunctionsDerivativeY(Ny, knots[kk].x(), knots[kk].y(), knots[kk].z());
1224 ShapeFunctionsDerivativeZ(Nz, knots[kk].x(), knots[kk].y(), knots[kk].z());
1225
1226 ddNz = m_ddT * Nz.transpose();
1227 d0d0Nz = m_d0d0T * Nz.transpose();
1228
1229 switch (kk) {
1230 case 0:
1231 case 1:
1232 case 2:
1233 case 3: {
1234 m_strainANS(kk) = 0.5 * ((Nz * ddNz)(0, 0) - (Nz * d0d0Nz)(0, 0));
1235 ChMatrixNM<double, 1, 3> tmpZ = Nz * m_d;
1236 for (int i = 0; i < 8; i++)
1237 for (int j = 0; j < 3; j++)
1238 m_strainANS_D(kk, i * 3 + j) = tmpZ(0, j) * Nz(0, i);
1239 break;
1240 }
1241 case 4:
1242 case 5: { // => yz
1243 m_strainANS(kk) = (Ny * ddNz)(0, 0) - (Ny * d0d0Nz)(0, 0);
1244 ChMatrixNM<double, 1, 3> tmpY = Ny * m_d;
1245 ChMatrixNM<double, 1, 3> tmpZ = Nz * m_d;
1246 for (int i = 0; i < 8; i++)
1247 for (int j = 0; j < 3; j++)
1248 m_strainANS_D(kk, i * 3 + j) = tmpY(0, j) * Nz(0, i) + tmpZ(0, j) * Ny(0, i);
1249 break;
1250 }
1251 case 6:
1252 case 7: { // => xz
1253 m_strainANS(kk) = (Nx * ddNz)(0, 0) - (Nx * d0d0Nz)(0, 0);
1254 ChMatrixNM<double, 1, 3> tmpX = Nx * m_d;
1255 ChMatrixNM<double, 1, 3> tmpZ = Nz * m_d;
1256 for (int i = 0; i < 8; i++)
1257 for (int j = 0; j < 3; j++)
1258 m_strainANS_D(kk, i * 3 + j) = tmpX(0, j) * Nz(0, i) + tmpZ(0, j) * Nx(0, i);
1259 break;
1260 }
1261 }
1262 }
1263 }
1264
1265 // -----------------------------------------------------------------------------
1266 // Interface to ChElementShell base class
1267 // -----------------------------------------------------------------------------
EvaluateDeflection(double & def)1268 void ChElementShellANCF_3423::EvaluateDeflection(double& def) {
1269 ChVector<> oldPos;
1270 ChVector<> newPos;
1271 ChVector<> defVec;
1272
1273 for (int i = 0; i < 4; i++) {
1274 oldPos.x() += this->m_d0(2 * i, 0);
1275 oldPos.y() += this->m_d0(2 * i, 1);
1276 oldPos.z() += this->m_d0(2 * i, 2);
1277 }
1278
1279 for (int i = 0; i < 4; i++) {
1280 newPos.x() += this->m_d(2 * i, 0);
1281 newPos.y() += this->m_d(2 * i, 1);
1282 newPos.z() += this->m_d(2 * i, 2);
1283 }
1284
1285 defVec = (newPos - oldPos) / 4;
1286 def = defVec.Length();
1287 }
1288
EvaluateSectionStrainStress(const ChVector<> & loc,int layer_id)1289 ChStrainStress3D ChElementShellANCF_3423::EvaluateSectionStrainStress(const ChVector<>& loc, int layer_id) {
1290 // Element shape function
1291 ShapeVector N;
1292 ShapeFunctions(N, loc.x(), loc.y(), loc.z());
1293
1294 // Determinant of position vector gradient matrix: Initial configuration
1295 ShapeVector Nx;
1296 ShapeVector Ny;
1297 ShapeVector Nz;
1298 ChMatrixNM<double, 1, 3> Nx_d0;
1299 ChMatrixNM<double, 1, 3> Ny_d0;
1300 ChMatrixNM<double, 1, 3> Nz_d0;
1301 double detJ0 = this->Calc_detJ0(loc.x(), loc.y(), loc.z(), Nx, Ny, Nz, Nx_d0, Ny_d0, Nz_d0);
1302
1303 // ANS shape function
1304 ChMatrixNM<double, 1, 4> S_ANS; // Shape function vector for Assumed Natural Strain
1305 ChMatrixNM<double, 6, 5> M; // Shape function vector for Enhanced Assumed Strain
1306 this->ShapeFunctionANSbilinearShell(S_ANS, loc.x(), loc.y());
1307 this->Basis_M(M, loc.x(), loc.y(), loc.z());
1308
1309 // Transformation : Orthogonal transformation (A and J)
1310 ChVector<double> G1xG2; // Cross product of first and second column of
1311 double G1dotG1; // Dot product of first column of position vector gradient
1312
1313 G1xG2.x() = Nx_d0(1) * Ny_d0(2) - Nx_d0(2) * Ny_d0(1);
1314 G1xG2.y() = Nx_d0(2) * Ny_d0(0) - Nx_d0(0) * Ny_d0(2);
1315 G1xG2.z() = Nx_d0(0) * Ny_d0(1) - Nx_d0(1) * Ny_d0(0);
1316 G1dotG1 = Nx_d0(0) * Nx_d0(0) + Nx_d0(1) * Nx_d0(1) + Nx_d0(2) * Nx_d0(2);
1317
1318 // Tangent Frame
1319 ChVector<double> A1;
1320 ChVector<double> A2;
1321 ChVector<double> A3;
1322 A1.x() = Nx_d0(0);
1323 A1.y() = Nx_d0(1);
1324 A1.z() = Nx_d0(2);
1325 A1 = A1 / sqrt(G1dotG1);
1326 A3 = G1xG2.GetNormalized();
1327 A2.Cross(A3, A1);
1328
1329 // Direction for orthotropic material
1330 ChVector<double> AA1;
1331 ChVector<double> AA2;
1332 ChVector<double> AA3;
1333 AA1 = A1;
1334 AA2 = A2;
1335 AA3 = A3;
1336
1337 /// Beta
1338 ChMatrixNM<double, 3, 3> j0;
1339 ChVector<double> j01;
1340 ChVector<double> j02;
1341 ChVector<double> j03;
1342 ChVectorN<double, 9> beta;
1343 // Calculates inverse of rd0 (j0) (position vector gradient: Initial Configuration)
1344 j0(0, 0) = Ny_d0(1) * Nz_d0(2) - Nz_d0(1) * Ny_d0(2);
1345 j0(0, 1) = Ny_d0(2) * Nz_d0(0) - Ny_d0(0) * Nz_d0(2);
1346 j0(0, 2) = Ny_d0(0) * Nz_d0(1) - Nz_d0(0) * Ny_d0(1);
1347 j0(1, 0) = Nz_d0(1) * Nx_d0(2) - Nx_d0(1) * Nz_d0(2);
1348 j0(1, 1) = Nz_d0(2) * Nx_d0(0) - Nx_d0(2) * Nz_d0(0);
1349 j0(1, 2) = Nz_d0(0) * Nx_d0(1) - Nz_d0(1) * Nx_d0(0);
1350 j0(2, 0) = Nx_d0(1) * Ny_d0(2) - Ny_d0(1) * Nx_d0(2);
1351 j0(2, 1) = Ny_d0(0) * Nx_d0(2) - Nx_d0(0) * Ny_d0(2);
1352 j0(2, 2) = Nx_d0(0) * Ny_d0(1) - Ny_d0(0) * Nx_d0(1);
1353 j0 /= detJ0;
1354
1355 j01[0] = j0(0, 0);
1356 j02[0] = j0(1, 0);
1357 j03[0] = j0(2, 0);
1358 j01[1] = j0(0, 1);
1359 j02[1] = j0(1, 1);
1360 j03[1] = j0(2, 1);
1361 j01[2] = j0(0, 2);
1362 j02[2] = j0(1, 2);
1363 j03[2] = j0(2, 2);
1364
1365 // Coefficients of contravariant transformation
1366 beta(0) = Vdot(AA1, j01);
1367 beta(1) = Vdot(AA2, j01);
1368 beta(2) = Vdot(AA3, j01);
1369 beta(3) = Vdot(AA1, j02);
1370 beta(4) = Vdot(AA2, j02);
1371 beta(5) = Vdot(AA3, j02);
1372 beta(6) = Vdot(AA1, j03);
1373 beta(7) = Vdot(AA2, j03);
1374 beta(8) = Vdot(AA3, j03);
1375
1376 // Transformation matrix, function of fiber angle
1377
1378 ChVectorN<double, 8> ddNx = m_ddT * Nx.transpose();
1379 ChVectorN<double, 8> ddNy = m_ddT * Ny.transpose();
1380
1381 ChVectorN<double, 8> d0d0Nx = m_d0d0T * Nx.transpose();
1382 ChVectorN<double, 8> d0d0Ny = m_d0d0T * Ny.transpose();
1383
1384 // Strain component
1385 ChVectorN<double, 6> strain_til;
1386 strain_til(0) = 0.5 * ((Nx * ddNx)(0, 0) - (Nx * d0d0Nx)(0, 0));
1387 strain_til(1) = 0.5 * ((Ny * ddNy)(0, 0) - (Ny * d0d0Ny)(0, 0));
1388 strain_til(2) = (Nx * ddNy)(0, 0) - (Nx * d0d0Ny)(0, 0);
1389 strain_til(3) = N(0) * this->m_strainANS(0) + N(2) * this->m_strainANS(1) + N(4) * this->m_strainANS(2) +
1390 N(6) * this->m_strainANS(3);
1391 strain_til(4) = S_ANS(0, 2) * this->m_strainANS(6) + S_ANS(0, 3) * this->m_strainANS(7);
1392 strain_til(5) = S_ANS(0, 0) * this->m_strainANS(4) + S_ANS(0, 1) * this->m_strainANS(5);
1393
1394 // For orthotropic material
1395 ChVectorN<double, 6> strain;
1396
1397 strain(0) = strain_til(0) * beta(0) * beta(0) + strain_til(1) * beta(3) * beta(3) +
1398 strain_til(2) * beta(0) * beta(3) + strain_til(3) * beta(6) * beta(6) +
1399 strain_til(4) * beta(0) * beta(6) + strain_til(5) * beta(3) * beta(6);
1400 strain(1) = strain_til(0) * beta(1) * beta(1) + strain_til(1) * beta(4) * beta(4) +
1401 strain_til(2) * beta(1) * beta(4) + strain_til(3) * beta(7) * beta(7) +
1402 strain_til(4) * beta(1) * beta(7) + strain_til(5) * beta(4) * beta(7);
1403 strain(2) = strain_til(0) * 2.0 * beta(0) * beta(1) + strain_til(1) * 2.0 * beta(3) * beta(4) +
1404 strain_til(2) * (beta(1) * beta(3) + beta(0) * beta(4)) + strain_til(3) * 2.0 * beta(6) * beta(7) +
1405 strain_til(4) * (beta(1) * beta(6) + beta(0) * beta(7)) +
1406 strain_til(5) * (beta(4) * beta(6) + beta(3) * beta(7));
1407 strain(3) = strain_til(0) * beta(2) * beta(2) + strain_til(1) * beta(5) * beta(5) +
1408 strain_til(2) * beta(2) * beta(5) + strain_til(3) * beta(8) * beta(8) +
1409 strain_til(4) * beta(2) * beta(8) + strain_til(5) * beta(5) * beta(8);
1410 strain(4) = strain_til(0) * 2.0 * beta(0) * beta(2) + strain_til(1) * 2.0 * beta(3) * beta(5) +
1411 strain_til(2) * (beta(2) * beta(3) + beta(0) * beta(5)) + strain_til(3) * 2.0 * beta(6) * beta(8) +
1412 strain_til(4) * (beta(2) * beta(6) + beta(0) * beta(8)) +
1413 strain_til(5) * (beta(5) * beta(6) + beta(3) * beta(8));
1414 strain(5) = strain_til(0) * 2.0 * beta(1) * beta(2) + strain_til(1) * 2.0 * beta(4) * beta(5) +
1415 strain_til(2) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3) * 2.0 * beta(7) * beta(8) +
1416 strain_til(4) * (beta(2) * beta(7) + beta(1) * beta(8)) +
1417 strain_til(5) * (beta(5) * beta(7) + beta(4) * beta(8));
1418
1419 const ChMatrixNM<double, 6, 6>& E_eps = GetLayer(layer_id).GetMaterial()->Get_E_eps();
1420 const ChVectorN<double, 6>& stress = E_eps * strain;
1421 ChStrainStress3D strainStressOut = {strain, stress};
1422
1423 return strainStressOut;
1424 }
EvaluateSectionDisplacement(const double u,const double v,ChVector<> & u_displ,ChVector<> & u_rotaz)1425 void ChElementShellANCF_3423::EvaluateSectionDisplacement(const double u,
1426 const double v,
1427 ChVector<>& u_displ,
1428 ChVector<>& u_rotaz) {
1429 // this is not a corotational element, so just do:
1430 EvaluateSectionPoint(u, v, u_displ);
1431 u_rotaz = VNULL; // no angles.. this is ANCF (or maybe return here the slope derivatives?)
1432 }
1433
EvaluateSectionFrame(const double u,const double v,ChVector<> & point,ChQuaternion<> & rot)1434 void ChElementShellANCF_3423::EvaluateSectionFrame(const double u,
1435 const double v,
1436 ChVector<>& point,
1437 ChQuaternion<>& rot) {
1438 // this is not a corotational element, so just do:
1439 EvaluateSectionPoint(u, v, point);
1440
1441 MatrixNx3 e_bar;
1442 CalcCoordMatrix(e_bar);
1443
1444 ShapeVector Nx;
1445 ShapeVector Ny;
1446 ShapeFunctionsDerivativeX(Nx, u, v, 0);
1447 ShapeFunctionsDerivativeY(Ny, u, v, 0);
1448
1449 // Since ANCF does not use rotations, calculate an approximate
1450 // rotation based off the position vector gradients
1451 ChVector<double> MidsurfaceX = e_bar.transpose() * Nx.transpose();
1452 ChVector<double> MidsurfaceY = e_bar.transpose() * Ny.transpose();
1453
1454 // Since the position vector gradients are not in general orthogonal,
1455 // set the Dx direction tangent to the shell xi axis and
1456 // compute the Dy and Dz directions by using a
1457 // Gram-Schmidt orthonormalization, guided by the shell eta axis
1458 ChMatrix33<> msect;
1459 msect.Set_A_Xdir(MidsurfaceX, MidsurfaceY);
1460
1461 rot = msect.Get_A_quaternion();
1462 }
1463
EvaluateSectionPoint(const double u,const double v,ChVector<> & point)1464 void ChElementShellANCF_3423::EvaluateSectionPoint(const double u, const double v, ChVector<>& point) {
1465 ChVector<> u_displ;
1466
1467 double x = u; // because ShapeFunctions() works in -1..1 range
1468 double y = v; // because ShapeFunctions() works in -1..1 range
1469 double z = 0;
1470 ShapeVector N;
1471 ShapeFunctions(N, x, y, z);
1472
1473 const ChVector<>& pA = m_nodes[0]->GetPos();
1474 const ChVector<>& pB = m_nodes[1]->GetPos();
1475 const ChVector<>& pC = m_nodes[2]->GetPos();
1476 const ChVector<>& pD = m_nodes[3]->GetPos();
1477
1478 point.x() = N(0) * pA.x() + N(2) * pB.x() + N(4) * pC.x() + N(6) * pD.x();
1479 point.y() = N(0) * pA.y() + N(2) * pB.y() + N(4) * pC.y() + N(6) * pD.y();
1480 point.z() = N(0) * pA.z() + N(2) * pB.z() + N(4) * pC.z() + N(6) * pD.z();
1481 }
1482
1483 // -----------------------------------------------------------------------------
1484 // Functions for ChLoadable interface
1485 // -----------------------------------------------------------------------------
1486
1487 // Gets all the DOFs packed in a single vector (position part).
LoadableGetStateBlock_x(int block_offset,ChState & mD)1488 void ChElementShellANCF_3423::LoadableGetStateBlock_x(int block_offset, ChState& mD) {
1489 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos().eigen();
1490 mD.segment(block_offset + 3, 3) = m_nodes[0]->GetD().eigen();
1491 mD.segment(block_offset + 6, 3) = m_nodes[1]->GetPos().eigen();
1492 mD.segment(block_offset + 9, 3) = m_nodes[1]->GetD().eigen();
1493 mD.segment(block_offset + 12, 3) = m_nodes[2]->GetPos().eigen();
1494 mD.segment(block_offset + 15, 3) = m_nodes[2]->GetD().eigen();
1495 mD.segment(block_offset + 18, 3) = m_nodes[3]->GetPos().eigen();
1496 mD.segment(block_offset + 21, 3) = m_nodes[3]->GetD().eigen();
1497 }
1498
1499 // Gets all the DOFs packed in a single vector (velocity part).
LoadableGetStateBlock_w(int block_offset,ChStateDelta & mD)1500 void ChElementShellANCF_3423::LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) {
1501 mD.segment(block_offset + 0, 3) = m_nodes[0]->GetPos_dt().eigen();
1502 mD.segment(block_offset + 3, 3) = m_nodes[0]->GetD_dt().eigen();
1503 mD.segment(block_offset + 6, 3) = m_nodes[1]->GetPos_dt().eigen();
1504 mD.segment(block_offset + 9, 3) = m_nodes[1]->GetD_dt().eigen();
1505 mD.segment(block_offset + 12, 3) = m_nodes[2]->GetPos_dt().eigen();
1506 mD.segment(block_offset + 15, 3) = m_nodes[2]->GetD_dt().eigen();
1507 mD.segment(block_offset + 18, 3) = m_nodes[3]->GetPos_dt().eigen();
1508 mD.segment(block_offset + 21, 3) = m_nodes[3]->GetD_dt().eigen();
1509 }
1510
LoadableStateIncrement(const unsigned int off_x,ChState & x_new,const ChState & x,const unsigned int off_v,const ChStateDelta & Dv)1511 void ChElementShellANCF_3423::LoadableStateIncrement(const unsigned int off_x,
1512 ChState& x_new,
1513 const ChState& x,
1514 const unsigned int off_v,
1515 const ChStateDelta& Dv) {
1516 for (int i = 0; i < 4; i++) {
1517 this->m_nodes[i]->NodeIntStateIncrement(off_x + 6 * i, x_new, x, off_v + 6 * i, Dv);
1518 }
1519 }
1520
EvaluateSectionVelNorm(double U,double V,ChVector<> & Result)1521 void ChElementShellANCF_3423::EvaluateSectionVelNorm(double U, double V, ChVector<>& Result) {
1522 ShapeVector N;
1523 ShapeFunctions(N, U, V, 0);
1524 for (unsigned int ii = 0; ii < 4; ii++) {
1525 Result += N(ii * 2) * m_nodes[ii]->GetPos_dt();
1526 Result += N(ii * 2 + 1) * m_nodes[ii]->GetPos_dt();
1527 }
1528 }
1529
1530 // Get the pointers to the contained ChVariables, appending to the mvars vector.
LoadableGetVariables(std::vector<ChVariables * > & mvars)1531 void ChElementShellANCF_3423::LoadableGetVariables(std::vector<ChVariables*>& mvars) {
1532 for (int i = 0; i < m_nodes.size(); ++i) {
1533 mvars.push_back(&m_nodes[i]->Variables());
1534 mvars.push_back(&m_nodes[i]->Variables_D());
1535 }
1536 }
1537
1538 // Evaluate N'*F , where N is the shape function evaluated at (U,V) coordinates of the surface.
1539 // This calculation takes a slightly different form for ANCF elements
1540 // For this ANCF element, only the first 6 entries in F are used in the calculation. The first three entries is
1541 // the applied force in global coordinates and the second 3 entries is the applied moment in global space.
ComputeNF(const double U,const double V,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1542 void ChElementShellANCF_3423::ComputeNF(
1543 const double U, // parametric coordinate in surface
1544 const double V, // parametric coordinate in surface
1545 ChVectorDynamic<>& Qi, // Return result of Q = N'*F here
1546 double& detJ, // Return det[J] here
1547 const ChVectorDynamic<>& F, // Input F vector, size is =n. field coords.
1548 ChVectorDynamic<>* state_x, // if != 0, update state (pos. part) to this, then evaluate Q
1549 ChVectorDynamic<>* state_w // if != 0, update state (speed part) to this, then evaluate Q
1550 ) {
1551 ShapeVector N;
1552 ShapeFunctions(N, U, V, 0);
1553
1554 Eigen::Map<MatrixNx3> QiCompact(Qi.data(), NSF, 3);
1555 QiCompact = N.transpose() * F.segment(0, 3).transpose();
1556
1557 // Compute the generalized force vector for the applied moment component
1558 // See: Antonio M Recuero, Javier F Aceituno, Jose L Escalona, and Ahmed A Shabana.
1559 // A nonlinear approach for modeling rail flexibility using the absolute nodal coordinate
1560 // formulation. Nonlinear Dynamics, 83(1-2):463-481, 2016.
1561 ShapeVector Nx;
1562 ShapeVector Ny;
1563 ShapeVector Nz;
1564 ShapeFunctionsDerivativeX(Nx, U, V, 0);
1565 ShapeFunctionsDerivativeY(Ny, U, V, 0);
1566 ShapeFunctionsDerivativeZ(Nz, U, V, 0);
1567
1568 // Change the shape function derivatives with respect to the element local "x", "y", or "z" coordinates to be with
1569 // respect to the normalized coordinates "xi", "eta", and "zeta" when loading the results into the Sxi_D matrix
1570 MatrixNx3 Sxi_D;
1571 Sxi_D.col(0) = Nx * m_lenX / 2;
1572 Sxi_D.col(1) = Ny * m_lenY / 2;
1573 Sxi_D.col(2) = Nz * m_thickness / 2;
1574
1575 MatrixNx3 e_bar;
1576 CalcCoordMatrix(e_bar);
1577
1578 // Calculate the element Jacobian between the current configuration and the normalized configuration
1579 ChMatrix33<double> J_Cxi = e_bar.transpose() * Sxi_D;
1580 ChMatrix33<double> J_Cxi_Inv = J_Cxi.inverse();
1581
1582 // Compute the unique pieces that make up the moment projection matrix "G"
1583 VectorN G_A = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 0) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 0) +
1584 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 0);
1585 VectorN G_B = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 1) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 1) +
1586 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 1);
1587 VectorN G_C = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 2) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 2) +
1588 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 2);
1589
1590 ChVectorN<double, 3> M_scaled = 0.5 * F.segment(3, 3);
1591
1592 // Compute G'M without actually forming the complete matrix "G" (since it has a sparsity pattern to it)
1593 for (unsigned int i = 0; i < NSF; i++) {
1594 Qi(3 * i + 0) += M_scaled(1) * G_C(i) - M_scaled(2) * G_B(i);
1595 Qi(3 * i + 1) += M_scaled(2) * G_A(i) - M_scaled(0) * G_C(i);
1596 Qi(3 * i + 2) += M_scaled(0) * G_B(i) - M_scaled(1) * G_A(i);
1597 }
1598
1599 // Compute the element Jacobian between the current configuration and the normalized configuration
1600 // This is different than the element Jacobian between the reference configuration and the normalized
1601 // configuration used in the internal force calculations. For this calculation, this is the ratio between the
1602 // actual differential area and the normalized differential area. The cross product is
1603 // used to calculate this area ratio for potential use in Gauss-Quadrature or similar numeric integration.
1604 detJ = J_Cxi.col(0).cross(J_Cxi.col(1)).norm();
1605 }
1606
1607 // Evaluate N'*F , where N is the shape function evaluated at (U,V,W) coordinates of the surface.
1608 // This calculation takes a slightly different form for ANCF elements
1609 // For this ANCF element, only the first 6 entries in F are used in the calculation. The first three entries is
1610 // the applied force in global coordinates and the second 3 entries is the applied moment in global space.
ComputeNF(const double U,const double V,const double W,ChVectorDynamic<> & Qi,double & detJ,const ChVectorDynamic<> & F,ChVectorDynamic<> * state_x,ChVectorDynamic<> * state_w)1611 void ChElementShellANCF_3423::ComputeNF(
1612 const double U, // parametric coordinate in volume
1613 const double V, // parametric coordinate in volume
1614 const double W, // parametric coordinate in volume
1615 ChVectorDynamic<>& Qi, // Return result of N'*F here, maybe with offset block_offset
1616 double& detJ, // Return det[J] here
1617 const ChVectorDynamic<>& F, // Input F vector, size is = n.field coords.
1618 ChVectorDynamic<>* state_x, // if != 0, update state (pos. part) to this, then evaluate Q
1619 ChVectorDynamic<>* state_w // if != 0, update state (speed part) to this, then evaluate Q
1620 ) {
1621 ShapeVector N;
1622 ShapeFunctions(N, U, V, W);
1623
1624 Eigen::Map<MatrixNx3> QiCompact(Qi.data(), NSF, 3);
1625 QiCompact = N.transpose() * F.segment(0, 3).transpose();
1626
1627 // Compute the generalized force vector for the applied moment component
1628 // See: Antonio M Recuero, Javier F Aceituno, Jose L Escalona, and Ahmed A Shabana.
1629 // A nonlinear approach for modeling rail flexibility using the absolute nodal coordinate
1630 // formulation. Nonlinear Dynamics, 83(1-2):463-481, 2016.
1631 ShapeVector Nx;
1632 ShapeVector Ny;
1633 ShapeVector Nz;
1634 ShapeFunctionsDerivativeX(Nx, U, V, W);
1635 ShapeFunctionsDerivativeY(Ny, U, V, W);
1636 ShapeFunctionsDerivativeZ(Nz, U, V, W);
1637
1638 // Change the shape function derivatives with respect to the element local "x", "y", or "z" coordinates to be with
1639 // respect to the normalized coordinates "xi", "eta", and "zeta" when loading the results into the Sxi_D matrix
1640 MatrixNx3 Sxi_D;
1641 Sxi_D.col(0) = Nx * m_lenX / 2;
1642 Sxi_D.col(1) = Ny * m_lenY / 2;
1643 Sxi_D.col(2) = Nz * m_thickness / 2;
1644
1645 MatrixNx3 e_bar;
1646 CalcCoordMatrix(e_bar);
1647
1648 // Calculate the element Jacobian between the current configuration and the normalized configuration
1649 ChMatrix33<double> J_Cxi = e_bar.transpose() * Sxi_D;
1650 ChMatrix33<double> J_Cxi_Inv = J_Cxi.inverse();
1651
1652 // Compute the unique pieces that make up the moment projection matrix "G"
1653 VectorN G_A = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 0) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 0) +
1654 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 0);
1655 VectorN G_B = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 1) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 1) +
1656 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 1);
1657 VectorN G_C = Sxi_D.col(0).transpose() * J_Cxi_Inv(0, 2) + Sxi_D.col(1).transpose() * J_Cxi_Inv(1, 2) +
1658 Sxi_D.col(2).transpose() * J_Cxi_Inv(2, 2);
1659
1660 ChVectorN<double, 3> M_scaled = 0.5 * F.segment(3, 3);
1661
1662 // Compute G'M without actually forming the complete matrix "G" (since it has a sparsity pattern to it)
1663 for (unsigned int i = 0; i < NSF; i++) {
1664 Qi(3 * i + 0) += M_scaled(1) * G_C(i) - M_scaled(2) * G_B(i);
1665 Qi(3 * i + 1) += M_scaled(2) * G_A(i) - M_scaled(0) * G_C(i);
1666 Qi(3 * i + 2) += M_scaled(0) * G_B(i) - M_scaled(1) * G_A(i);
1667 }
1668
1669 // Compute the element Jacobian between the current configuration and the normalized configuration
1670 // This is different than the element Jacobian between the reference configuration and the normalized
1671 // configuration used in the internal force calculations. For this calculation, this is the ratio between the
1672 // actual differential volume and the normalized differential volume. The determinate of the element Jacobian is
1673 // used to calculate this volume ratio for potential use in Gauss-Quadrature or similar numeric integration.
1674 detJ = J_Cxi.determinant();
1675 }
1676
1677 // -----------------------------------------------------------------------------
1678 //
1679 // -----------------------------------------------------------------------------
1680
1681 // Calculate average element density (needed for ChLoaderVolumeGravity).
GetDensity()1682 double ChElementShellANCF_3423::GetDensity() {
1683 double tot_density = 0;
1684 for (size_t kl = 0; kl < m_numLayers; kl++) {
1685 double rho = m_layers[kl].GetMaterial()->Get_rho();
1686 double layerthick = m_layers[kl].Get_thickness();
1687 tot_density += rho * layerthick;
1688 }
1689 return tot_density / m_thickness;
1690 }
1691
1692 // Calculate normal to the surface at (U,V) coordinates.
ComputeNormal(const double U,const double V)1693 ChVector<> ChElementShellANCF_3423::ComputeNormal(const double U, const double V) {
1694 ChMatrixNM<double, 8, 3> mD;
1695 ChMatrixNM<double, 1, 8> Nx;
1696 ChMatrixNM<double, 1, 8> Ny;
1697 ChMatrixNM<double, 1, 8> Nz;
1698
1699 ShapeFunctionsDerivativeX(Nx, U, V, 0);
1700 ShapeFunctionsDerivativeY(Ny, U, V, 0);
1701 ShapeFunctionsDerivativeZ(Nz, U, V, 0);
1702
1703 CalcCoordMatrix(mD);
1704
1705 ChMatrixNM<double, 1, 3> Nx_d = Nx * mD;
1706 ChMatrixNM<double, 1, 3> Ny_d = Ny * mD;
1707 ChMatrixNM<double, 1, 3> Nz_d = Nz * mD;
1708
1709 ChMatrixNM<double, 3, 3> rd;
1710 rd(0, 0) = Nx_d(0, 0);
1711 rd(1, 0) = Nx_d(0, 1);
1712 rd(2, 0) = Nx_d(0, 2);
1713 rd(0, 1) = Ny_d(0, 0);
1714 rd(1, 1) = Ny_d(0, 1);
1715 rd(2, 1) = Ny_d(0, 2);
1716 rd(0, 2) = Nz_d(0, 0);
1717 rd(1, 2) = Nz_d(0, 1);
1718 rd(2, 2) = Nz_d(0, 2);
1719
1720 ChVector<> G1xG2;
1721 G1xG2[0] = rd(1, 0) * rd(2, 1) - rd(2, 0) * rd(1, 1);
1722 G1xG2[1] = rd(2, 0) * rd(0, 1) - rd(0, 0) * rd(2, 1);
1723 G1xG2[2] = rd(0, 0) * rd(1, 1) - rd(1, 0) * rd(0, 1);
1724
1725 double G1xG2nrm = sqrt(G1xG2[0] * G1xG2[0] + G1xG2[1] * G1xG2[1] + G1xG2[2] * G1xG2[2]);
1726 return G1xG2 / G1xG2nrm;
1727 }
1728
1729 // ============================================================================
1730 // Implementation of ChElementShellANCF_3423::Layer methods
1731 // ============================================================================
1732
1733 // Private constructor (a layer can be created only by adding it to an element)
Layer(ChElementShellANCF_3423 * element,double thickness,double theta,std::shared_ptr<ChMaterialShellANCF> material)1734 ChElementShellANCF_3423::Layer::Layer(ChElementShellANCF_3423* element,
1735 double thickness,
1736 double theta,
1737 std::shared_ptr<ChMaterialShellANCF> material)
1738 : m_element(element), m_thickness(thickness), m_theta(theta), m_material(material) {}
1739
1740 // Initial setup for this layer: calculate T0 and detJ0 at the element center.
SetupInitial()1741 void ChElementShellANCF_3423::Layer::SetupInitial() {
1742 // Evaluate shape functions at element center
1743 ChMatrixNM<double, 1, 8> Nx;
1744 ChMatrixNM<double, 1, 8> Ny;
1745 ChMatrixNM<double, 1, 8> Nz;
1746 m_element->ShapeFunctionsDerivativeX(Nx, 0, 0, 0);
1747 m_element->ShapeFunctionsDerivativeY(Ny, 0, 0, 0);
1748 m_element->ShapeFunctionsDerivativeZ(Nz, 0, 0, 0);
1749
1750 ChMatrixNM<double, 1, 3> Nx_d0 = Nx * m_element->m_d0;
1751 ChMatrixNM<double, 1, 3> Ny_d0 = Ny * m_element->m_d0;
1752 ChMatrixNM<double, 1, 3> Nz_d0 = Nz * m_element->m_d0;
1753
1754 // Determinant of position vector gradient matrix: Initial configuration
1755 m_detJ0C = Nx_d0(0, 0) * Ny_d0(0, 1) * Nz_d0(0, 2) + Ny_d0(0, 0) * Nz_d0(0, 1) * Nx_d0(0, 2) +
1756 Nz_d0(0, 0) * Nx_d0(0, 1) * Ny_d0(0, 2) - Nx_d0(0, 2) * Ny_d0(0, 1) * Nz_d0(0, 0) -
1757 Ny_d0(0, 2) * Nz_d0(0, 1) * Nx_d0(0, 0) - Nz_d0(0, 2) * Nx_d0(0, 1) * Ny_d0(0, 0);
1758
1759 //// Transformation : Orthogonal transformation (A and J) ////
1760 ChVector<double> G1xG2; // Cross product of first and second column of
1761 double G1dotG1; // Dot product of first column of position vector gradient
1762
1763 G1xG2.x() = Nx_d0(1) * Ny_d0(2) - Nx_d0(2) * Ny_d0(1);
1764 G1xG2.y() = Nx_d0(2) * Ny_d0(0) - Nx_d0(0) * Ny_d0(2);
1765 G1xG2.z() = Nx_d0(0) * Ny_d0(1) - Nx_d0(1) * Ny_d0(0);
1766 G1dotG1 = Nx_d0(0) * Nx_d0(0) + Nx_d0(1) * Nx_d0(1) + Nx_d0(2) * Nx_d0(2);
1767
1768 // Tangent Frame
1769 ChVector<double> A1;
1770 ChVector<double> A2;
1771 ChVector<double> A3;
1772 A1.x() = Nx_d0(0);
1773 A1.y() = Nx_d0(1);
1774 A1.z() = Nx_d0(2);
1775 A1 = A1 / sqrt(G1dotG1);
1776 A3 = G1xG2.GetNormalized();
1777 A2.Cross(A3, A1);
1778
1779 ChVector<double> AA1;
1780 ChVector<double> AA2;
1781 ChVector<double> AA3;
1782 AA1 = A1 * cos(m_theta) + A2 * sin(m_theta);
1783 AA2 = -A1 * sin(m_theta) + A2 * cos(m_theta);
1784 AA3 = A3;
1785
1786 ////Beta
1787 ChMatrixNM<double, 3, 3> j0;
1788 ChVector<double> j01;
1789 ChVector<double> j02;
1790 ChVector<double> j03;
1791 ChVectorN<double, 9> beta;
1792
1793 j0(0, 0) = Ny_d0(1) * Nz_d0(2) - Nz_d0(1) * Ny_d0(2);
1794 j0(0, 1) = Ny_d0(2) * Nz_d0(0) - Ny_d0(0) * Nz_d0(2);
1795 j0(0, 2) = Ny_d0(0) * Nz_d0(1) - Nz_d0(0) * Ny_d0(1);
1796 j0(1, 0) = Nz_d0(1) * Nx_d0(2) - Nx_d0(1) * Nz_d0(2);
1797 j0(1, 1) = Nz_d0(2) * Nx_d0(0) - Nx_d0(2) * Nz_d0(0);
1798 j0(1, 2) = Nz_d0(0) * Nx_d0(1) - Nz_d0(1) * Nx_d0(0);
1799 j0(2, 0) = Nx_d0(1) * Ny_d0(2) - Ny_d0(1) * Nx_d0(2);
1800 j0(2, 1) = Ny_d0(0) * Nx_d0(2) - Nx_d0(0) * Ny_d0(2);
1801 j0(2, 2) = Nx_d0(0) * Ny_d0(1) - Ny_d0(0) * Nx_d0(1);
1802 j0 /= m_detJ0C;
1803
1804 j01[0] = j0(0, 0);
1805 j02[0] = j0(1, 0);
1806 j03[0] = j0(2, 0);
1807 j01[1] = j0(0, 1);
1808 j02[1] = j0(1, 1);
1809 j03[1] = j0(2, 1);
1810 j01[2] = j0(0, 2);
1811 j02[2] = j0(1, 2);
1812 j03[2] = j0(2, 2);
1813
1814 beta(0) = Vdot(AA1, j01);
1815 beta(1) = Vdot(AA2, j01);
1816 beta(2) = Vdot(AA3, j01);
1817 beta(3) = Vdot(AA1, j02);
1818 beta(4) = Vdot(AA2, j02);
1819 beta(5) = Vdot(AA3, j02);
1820 beta(6) = Vdot(AA1, j03);
1821 beta(7) = Vdot(AA2, j03);
1822 beta(8) = Vdot(AA3, j03);
1823
1824 // Calculate T0: transformation matrix, function of fiber angle (see Yamashita et al, 2015, JCND)
1825 m_T0(0, 0) = pow(beta(0), 2);
1826 m_T0(1, 0) = pow(beta(1), 2);
1827 m_T0(2, 0) = 2.0 * beta(0) * beta(1);
1828 m_T0(3, 0) = pow(beta(2), 2);
1829 m_T0(4, 0) = 2.0 * beta(0) * beta(2);
1830 m_T0(5, 0) = 2.0 * beta(1) * beta(2);
1831
1832 m_T0(0, 1) = pow(beta(3), 2);
1833 m_T0(1, 1) = pow(beta(4), 2);
1834 m_T0(2, 1) = 2.0 * beta(3) * beta(4);
1835 m_T0(3, 1) = pow(beta(5), 2);
1836 m_T0(4, 1) = 2.0 * beta(3) * beta(5);
1837 m_T0(5, 1) = 2.0 * beta(4) * beta(5);
1838
1839 m_T0(0, 2) = beta(0) * beta(3);
1840 m_T0(1, 2) = beta(1) * beta(4);
1841 m_T0(2, 2) = beta(0) * beta(4) + beta(1) * beta(3);
1842 m_T0(3, 2) = beta(2) * beta(5);
1843 m_T0(4, 2) = beta(0) * beta(5) + beta(2) * beta(3);
1844 m_T0(5, 2) = beta(2) * beta(4) + beta(1) * beta(5);
1845
1846 m_T0(0, 3) = pow(beta(6), 2);
1847 m_T0(1, 3) = pow(beta(7), 2);
1848 m_T0(2, 3) = 2.0 * beta(6) * beta(7);
1849 m_T0(3, 3) = pow(beta(8), 2);
1850 m_T0(4, 3) = 2.0 * beta(6) * beta(8);
1851 m_T0(5, 3) = 2.0 * beta(7) * beta(8);
1852
1853 m_T0(0, 4) = beta(0) * beta(6);
1854 m_T0(1, 4) = beta(1) * beta(7);
1855 m_T0(2, 4) = beta(0) * beta(7) + beta(6) * beta(1);
1856 m_T0(3, 4) = beta(2) * beta(8);
1857 m_T0(4, 4) = beta(0) * beta(8) + beta(2) * beta(6);
1858 m_T0(5, 4) = beta(1) * beta(8) + beta(2) * beta(7);
1859
1860 m_T0(0, 5) = beta(3) * beta(6);
1861 m_T0(1, 5) = beta(4) * beta(7);
1862 m_T0(2, 5) = beta(3) * beta(7) + beta(4) * beta(6);
1863 m_T0(3, 5) = beta(5) * beta(8);
1864 m_T0(4, 5) = beta(3) * beta(8) + beta(6) * beta(5);
1865 m_T0(5, 5) = beta(4) * beta(8) + beta(5) * beta(7);
1866 }
1867
1868 } // end of namespace fea
1869 } // end of namespace chrono
1870