1 /** @file gsMaterialMatrix.hpp
2
3 @brief Provides hyperelastic material matrices
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s):
12 H.M. Verhelst (2019-..., TU Delft)
13 A. Mantzaflaris (2019-..., Inria)
14 */
15
16 /*
17 To Do [updated 16-06-2020]:
18 - Make beta (compressible materials) and material parameters universal for all integration points over the thickness. So get them out of the _dPsi functions etc and move them into the integration loops as global variables.
19
20 */
21
22
23
24 #pragma once
25
26 #include <gsKLShell/gsMaterialMatrix.h>
27
28 namespace gismo
29 {
30
31 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
gsMaterialMatrix(const gsFunctionSet<T> & mp,const gsFunction<T> & thickness,const std::vector<gsFunction<T> * > & pars)32 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::gsMaterialMatrix(
33 const gsFunctionSet<T> & mp,
34 const gsFunction<T> & thickness,
35 const std::vector<gsFunction<T>*> &pars
36 )
37 :
38 Base(mp),
39 m_thickness(&thickness),
40 m_pars(pars)
41 {
42 _initialize();
43 }
44
45 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
gsMaterialMatrix(const gsFunctionSet<T> & mp,const gsFunction<T> & thickness,const std::vector<gsFunction<T> * > & pars,const gsFunction<T> & density)46 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::gsMaterialMatrix(
47 const gsFunctionSet<T> & mp,
48 const gsFunction<T> & thickness,
49 const std::vector<gsFunction<T>*> &pars,
50 const gsFunction<T> & density
51 )
52 :
53 Base(mp),
54 m_thickness(&thickness),
55 m_pars(pars),
56 m_density(&density)
57 {
58 _initialize();
59 }
60
61 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
info() const62 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::info() const
63 {
64 gsInfo <<"---------------------------------------------------------------------\n"
65 <<"---------------------Hyperelastic Material Info----------------------\n"
66 <<"---------------------------------------------------------------------\n\n";
67
68 gsInfo <<"Material model: \t";
69 if (comp)
70 gsInfo<<"Compressible ";
71 else
72 gsInfo<<"Incompressible ";
73
74 if (mat==Material::SvK)
75 gsInfo<<"Saint-Venant Kirchhoff";
76 else if (mat==Material::NH)
77 gsInfo<<"Neo-Hookean";
78 else if (mat==Material::MR)
79 gsInfo<<"Mooney-Rivlin";
80 else if (mat==Material::OG)
81 gsInfo<<"Ogden";
82 else if (mat==Material::NH_ext)
83 gsInfo<<"Neo-Hookean Extended";
84 else
85 gsWarn<<"Not specified";
86 gsInfo<<"\n";
87
88 gsInfo <<"Implementation: \t";
89 if (imp==Implementation::Analytical)
90 gsInfo<<"Analytical";
91 else if (imp==Implementation::Generalized)
92 gsInfo<<"Generalized";
93 else if (imp==Implementation::Spectral)
94 gsInfo<<"Spectral";
95 else
96 gsWarn<<"Not specified";
97 gsInfo<<" implementation\n";
98
99 gsInfo <<"---------------------------------------------------------------------\n\n";
100
101 }
102
103 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_defaultOptions()104 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_defaultOptions()
105 {
106 m_options.addInt("NumGauss","Number of Gaussian points through thickness",4);
107 }
108
109 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_initialize()110 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_initialize()
111 {
112 // Set default options
113 this->_defaultOptions();
114
115 // set flags
116 m_map.mine().flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
117
118 // Initialize some parameters
119 m_numPars = m_pars.size();
120 }
121
122
123 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_computePoints(const index_t patch,const gsMatrix<T> & u) const124 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_computePoints(const index_t patch, const gsMatrix<T> & u) const
125 {
126 gsMatrix<T> tmp;
127
128 this->_computeMetricUndeformed(patch,u);
129
130 if (Base::m_defpatches->nPieces()!=0)
131 this->_computeMetricDeformed(patch,u);
132
133 m_thickness->eval_into(m_map.mine().values[0], m_Tmat);
134
135 m_parmat.resize(m_numPars,m_map.mine().values[0].cols());
136 m_parmat.setZero();
137
138 for (size_t v=0; v!=m_pars.size(); v++)
139 {
140 m_pars[v]->eval_into(m_map.mine().values[0], tmp);
141 m_parmat.row(v) = tmp;
142 }
143
144 m_parvals.resize(m_numPars);
145
146 _computePoints_impl<mat>(u);
147 }
148
149 // CHECK ONLY
150 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
151 template <enum Material _mat>
152 typename std::enable_if<_mat==Material::OG, void>::type
_computePoints_impl(const gsMatrix<T> & u) const153 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_computePoints_impl(const gsMatrix<T>& u) const
154 {
155 // T prod, sum, mu;
156 for (index_t c=0; c!=m_parmat.cols(); c++)
157 {
158 for (index_t r=0; r!=m_parmat.rows(); r++)
159 m_parvals.at(r) = m_parmat(r,c);
160
161 // mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
162 GISMO_ENSURE((m_numPars-2 )% 2 ==0, "Ogden material models must have an even number of parameters (tuples of alpha_i and mu_i). m_numPars = "<< m_numPars);
163
164 /// THIS CHECK IS NOT NECESSARY
165 // int n = (m_numPars-2)/2;
166 // sum = 0.0;
167 // for (index_t k=0; k!=n; k++)
168 // {
169 // prod = m_parvals.at(2*(k+1))*m_parvals.at(2*(k+1)+1);
170 // GISMO_ENSURE(prod > 0.0,"Product of coefficients must be positive for all indices");
171 // sum += prod;
172 // }
173 // GISMO_ENSURE((sum-2.*mu)/sum<1e-10,"Sum of products must be equal to 2*mu! sum = "<<sum<<"; 2*mu = "<<2.*mu);
174 }
175 }
176
177 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
178 template <enum Material _mat>
179 typename std::enable_if<_mat!=Material::OG, void>::type
_computePoints_impl(const gsMatrix<T> & u) const180 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_computePoints_impl(const gsMatrix<T>& u) const
181 {
182
183 }
184
185 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
density_into(const index_t patch,const gsMatrix<T> & u,gsMatrix<T> & result) const186 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::density_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
187 {
188 m_map.mine().flags = NEED_VALUE;
189 m_map.mine().points = u;
190 static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(m_map);
191
192 result.resize(1, u.cols());
193 m_thickness->eval_into(m_map.mine().values[0], m_Tmat);
194 m_density->eval_into(m_map.mine().values[0], m_rhomat);
195 for (index_t i = 0; i != u.cols(); ++i) // points
196 {
197 result(0,i) = m_Tmat(0,i)*m_rhomat(0,i);
198 }
199
200 }
201
202 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
stretch_into(const index_t patch,const gsMatrix<T> & u,gsMatrix<T> & result) const203 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::stretch_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
204 {
205 m_map.mine().points = u;
206 static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(m_map);
207
208 this->_computePoints(patch,u);
209
210 _stretch_into_impl<comp>(u,result);
211 }
212
213 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
214 template <bool _comp>
215 typename std::enable_if<!_comp, void>::type
_stretch_into_impl(const gsMatrix<T> & u,gsMatrix<T> & result) const216 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_stretch_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
217 {
218 result.resize(3, u.cols());
219 std::pair<gsVector<T>,gsMatrix<T>> res;
220 for (index_t i=0; i!= u.cols(); i++)
221 {
222 this->_getMetric(i,0.0); // on point i, with height 0.0
223
224 gsMatrix<T> C(3,3);
225 C.setZero();
226 C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
227 C(2,2) = 1./m_J0_sq;
228
229 res = this->_evalStretch(C);
230 result.col(i) = res.first;
231 }
232 }
233
234 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
235 template <bool _comp>
236 typename std::enable_if<_comp, void>::type
_stretch_into_impl(const gsMatrix<T> & u,gsMatrix<T> & result) const237 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_stretch_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
238 {
239 result.resize(3, u.cols());
240 std::pair<gsVector<T>,gsMatrix<T>> res;
241 for (index_t i=0; i!= u.cols(); i++)
242 {
243 for (index_t v=0; v!=m_parmat.rows(); v++)
244 m_parvals.at(v) = m_parmat(v,i);
245
246 this->_getMetric(i,0.0); // on point i, with height 0.0
247
248 // Define objects
249 gsMatrix<T,3,3> c, cinv;
250 T S33, C3333, dc33;
251 // T S33_old;
252 S33 = 0.0;
253 dc33 = 0.0;
254 C3333 = 1.0;
255
256 index_t itmax = 100;
257 T tol = 1e-10;
258
259 // Initialize c
260 c.setZero();
261 c.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
262 c(2,2) = math::pow(m_J0_sq,-1.0); // c33
263 // c(2,2) = 1.0; // c33
264 cinv.block(0,0,2,2) = m_Gcon_def.block(0,0,2,2);
265 cinv(2,2) = 1.0/c(2,2);
266
267 m_J_sq = m_J0_sq * c(2,2);
268 S33 = _Sij(2,2,c,cinv);
269 // S33_old = (S33 == 0.0) ? 1.0 : S33;
270 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
271
272 dc33 = -2. * S33 / C3333;
273 for (index_t it = 0; it < itmax; it++)
274 {
275 c(2,2) += dc33;
276
277 GISMO_ENSURE(c(2,2)>= 0,"ERROR! c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
278 cinv(2,2) = 1.0/c(2,2);
279
280 m_J_sq = m_J0_sq * c(2,2) ;
281
282 S33 = _Sij(2,2,c,cinv);
283 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
284
285 dc33 = -2. * S33 / C3333;
286 if (abs(dc33) < tol)
287 {
288 res = this->_evalStretch(c);
289 result.col(i) = res.first;
290 break;
291 }
292 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, abs(dc33) = "<<abs(dc33)<<" and tolerance = "<<tol<<"\n");
293 }
294 }
295 }
296
297
298 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
stretchDir_into(const index_t patch,const gsMatrix<T> & u,gsMatrix<T> & result) const299 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::stretchDir_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
300 {
301 m_map.mine().points = u;
302 static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(m_map);
303
304 this->_computePoints(patch,u);
305
306 _stretchDir_into_impl<comp>(u,result);
307 }
308
309 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
310 template <bool _comp>
311 typename std::enable_if<!_comp, void>::type
_stretchDir_into_impl(const gsMatrix<T> & u,gsMatrix<T> & result) const312 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_stretchDir_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
313 {
314 result.resize(9, u.cols());
315 std::pair<gsVector<T>,gsMatrix<T>> res;
316 for (index_t i=0; i!= u.cols(); i++)
317 {
318 this->_getMetric(i,0.0); // on point i, with height 0.0
319
320 gsMatrix<T> C(3,3);
321 C.setZero();
322 C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
323 C(2,2) = 1./m_J0_sq;
324
325 res = this->_evalStretch(C);
326 result.col(i) = res.second.reshape(9,1);
327 }
328 }
329
330 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
331 template <bool _comp>
332 typename std::enable_if<_comp, void>::type
_stretchDir_into_impl(const gsMatrix<T> & u,gsMatrix<T> & result) const333 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_stretchDir_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
334 {
335 result.resize(9, u.cols());
336 std::pair<gsVector<T>,gsMatrix<T>> res;
337 for (index_t i=0; i!= u.cols(); i++)
338 {
339 for (index_t v=0; v!=m_parmat.rows(); v++)
340 m_parvals.at(v) = m_parmat(v,i);
341
342 this->_getMetric(i,0.0); // on point i, with height 0.0
343
344 // Define objects
345 gsMatrix<T,3,3> c, cinv;
346 T S33, C3333, dc33;
347 // T S33_old;
348 S33 = 0.0;
349 dc33 = 0.0;
350 C3333 = 1.0;
351
352 index_t itmax = 100;
353 T tol = 1e-10;
354
355 // Initialize c
356 c.setZero();
357 c.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
358 c(2,2) = math::pow(m_J0_sq,-1.0); // c33
359 // c(2,2) = 1.0; // c33
360 cinv.block(0,0,2,2) = m_Gcon_def.block(0,0,2,2);
361 cinv(2,2) = 1.0/c(2,2);
362
363 m_J_sq = m_J0_sq * c(2,2);
364 S33 = _Sij(2,2,c,cinv);
365 // S33_old = (S33 == 0.0) ? 1.0 : S33;
366 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
367
368 dc33 = -2. * S33 / C3333;
369 for (index_t it = 0; it < itmax; it++)
370 {
371 c(2,2) += dc33;
372
373 GISMO_ENSURE(c(2,2)>= 0,"ERROR! c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
374 cinv(2,2) = 1.0/c(2,2);
375
376 m_J_sq = m_J0_sq * c(2,2) ;
377
378 S33 = _Sij(2,2,c,cinv);
379 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
380
381 dc33 = -2. * S33 / C3333;
382 if (abs(dc33) < tol)
383 {
384 res = this->_evalStretch(c);
385 result.col(i) = res.second.reshape(9,1);
386 }
387 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, abs(dc33) = "<<abs(dc33)<<" and tolerance = "<<tol<<"\n");
388 }
389 }
390 }
391
392 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
thickness_into(const index_t patch,const gsMatrix<T> & u,gsMatrix<T> & result) const393 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::thickness_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
394 {
395 m_map.mine().flags = NEED_VALUE;
396 m_map.mine().points = u;
397 static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(m_map);
398 m_thickness->eval_into(m_map.mine().values[0], result);
399 }
400
401 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
transform_into(const index_t patch,const gsMatrix<T> & u,gsMatrix<T> & result) const402 void gsMaterialMatrix<dim,T,matId,comp,mat,imp>::transform_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
403 {
404 result.resize(9, u.cols());
405 gsMatrix<T> tmp, conbasis,sbasis;
406 this->stretchDir_into(patch,u,tmp);
407 for (index_t i=0; i!= u.cols(); i++)
408 {
409 this->_getMetric(i,0.0,true); // on point i, with height 0.0
410 sbasis = tmp.reshapeCol(i,3,3);
411 conbasis = m_gcon_ori;
412 result.col(i) = this->_transformation(conbasis,sbasis).reshape(9,1);
413 }
414 }
415
416 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
eval3D_matrix(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z,enum MaterialOutput out) const417 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::eval3D_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
418 {
419 return _eval3D_matrix_impl<mat,comp>(patch,u,z);
420 }
421
422 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
423 template <enum Material _mat, bool _comp>
424 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
_eval3D_matrix_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const425 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
426 {
427 gsMatrix<T> result = _eval_Incompressible_matrix(patch, u, z);
428 return result;
429 }
430
431 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
432 template <enum Material _mat, bool _comp>
433 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_matrix_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const434 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
435 {
436 gsMatrix<T> result = _eval_Incompressible_matrix(patch, u, z);
437 return result;
438 }
439
440 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
441 template <enum Material _mat, bool _comp>
442 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_matrix_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const443 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
444 {
445 gsMatrix<T> result = _eval_Compressible_matrix(patch, u, z);
446 return result;
447 }
448
449 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
eval3D_vector(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z,enum MaterialOutput out) const450 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::eval3D_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
451 {
452 return _eval3D_vector_impl<mat,comp>(patch,u,z);
453 }
454
455 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
456 template <enum Material _mat, bool _comp>
457 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
_eval3D_vector_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const458 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_vector_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
459 {
460 gsMatrix<T> result = _eval_Incompressible_vector(patch, u, z);
461 return result;
462 }
463
464 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
465 template <enum Material _mat, bool _comp>
466 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_vector_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const467 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_vector_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
468 {
469 gsMatrix<T> result = _eval_Incompressible_vector(patch, u, z);
470 return result;
471 }
472
473 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
474 template <enum Material _mat, bool _comp>
475 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_vector_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const476 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_vector_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
477 {
478 gsMatrix<T> result = _eval_Compressible_vector(patch, u, z);
479 return result;
480 }
481
482 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
eval3D_pstress(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z,enum MaterialOutput out) const483 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::eval3D_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
484 {
485 return _eval3D_pstress_impl<mat,comp>(patch,u,z);
486 }
487
488 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
489 template <enum Material _mat, bool _comp>
490 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
_eval3D_pstress_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const491 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_pstress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
492 {
493 gsMatrix<T> result = _eval_Incompressible_pstress(patch, u, z);
494 return result;
495 }
496
497 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
498 template <enum Material _mat, bool _comp>
499 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_pstress_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const500 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_pstress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
501 {
502 gsMatrix<T> result = _eval_Incompressible_pstress(patch, u, z);
503 return result;
504 }
505
506 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
507 template <enum Material _mat, bool _comp>
508 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
_eval3D_pstress_impl(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const509 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval3D_pstress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
510 {
511 gsMatrix<T> result = _eval_Compressible_pstress(patch, u, z);
512 return result;
513 }
514
515 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Incompressible_matrix(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const516 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Incompressible_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
517 {
518 // gsInfo<<"TO DO: evaluate moments using thickness";
519 // Input: u in-plane points
520 // z matrix with, per point, a column with z integration points
521 // Output: (n=u.cols(), m=z.rows())
522 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
523
524 this->_computePoints(patch,u);
525 gsMatrix<T> result(9, u.cols() * z.rows());
526 result.setZero();
527
528 for (index_t k=0; k!=u.cols(); k++)
529 {
530 // Evaluate material properties on the quadrature point
531 for (index_t v=0; v!=m_parmat.rows(); v++)
532 m_parvals.at(v) = m_parmat(v,k);
533
534 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
535 {
536 // this->computeMetric(i,z.at(j),true,true); // on point i, on height z(0,j)
537 this->_getMetric(k,z(j,k) * m_Tmat(0,k) ); // on point i, on height z(0,j)
538
539 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(j*u.cols()+k,3,3);
540 /*
541 C = C1111, C1122, C1112
542 symm, C2222, C2212
543 symm, symm, C1212
544 */
545 C(0,0) = _Cijkl(0,0,0,0); // C1111
546 C(1,1) = _Cijkl(1,1,1,1); // C2222
547 C(2,2) = _Cijkl(0,1,0,1); // C1212
548 C(1,0) = C(0,1) = _Cijkl(0,0,1,1); // C1122
549 C(2,0) = C(0,2) = _Cijkl(0,0,0,1); // C1112
550 C(2,1) = C(1,2) = _Cijkl(1,1,0,1); // C2212
551 }
552 }
553
554 return result;
555 }
556
557 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Incompressible_vector(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const558 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Incompressible_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
559 {
560 // gsInfo<<"TO DO: evaluate moments using thickness";
561 // Input: u in-plane points
562 // z matrix with, per point, a column with z integration points
563 // Output: (n=u.cols(), m=z.rows())
564 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
565
566 this->_computePoints(patch,u);
567 gsMatrix<T> result(3, u.cols() * z.rows());
568 result.setZero();
569
570 for (index_t k=0; k!=u.cols(); k++)
571 {
572 // Evaluate material properties on the quadrature point
573 for (index_t v=0; v!=m_parmat.rows(); v++)
574 m_parvals.at(v) = m_parmat(v,k);
575
576 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
577 {
578 this->_getMetric(k, z(j, k) * m_Tmat(0, k)); // on point i, on height z(0,j)
579
580 result(0, j * u.cols() + k) = _Sij(0, 0);
581 result(1, j * u.cols() + k) = _Sij(1, 1);
582 result(2, j * u.cols() + k) = _Sij(0, 1);
583 }
584 }
585
586 return result;
587 }
588
589 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Incompressible_pstress(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const590 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Incompressible_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
591 {
592 // gsInfo<<"TO DO: evaluate moments using thickness";
593 // Input: u in-plane points
594 // z matrix with, per point, a column with z integration points
595 // Output: (n=u.cols(), m=z.rows())
596 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
597
598 this->_computePoints(patch,u);
599 gsMatrix<T> result(3, u.cols() * z.rows());
600 result.setZero();
601
602 for (index_t k=0; k!=u.cols(); k++)
603 {
604 // Evaluate material properties on the quadrature point
605 for (index_t v=0; v!=m_parmat.rows(); v++)
606 m_parvals.at(v) = m_parmat(v,k);
607
608 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
609 {
610 this->_getMetric(k, z(j, k) * m_Tmat(0, k)); // on point i, on height z(0,j)
611
612 result(0, j * u.cols() + k) = _Sii(0);
613 result(1, j * u.cols() + k) = _Sii(1);
614 result(2, j * u.cols() + k) = 0;
615 }
616 }
617
618 return result;
619 }
620
621 /*
622 Available class members:
623 - m_parvals
624 - m_metric
625 - m_metric_def
626 - m_J0
627 - m_J
628 */
629 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
630 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl(const index_t i, const index_t j, const index_t k, const index_t l) const
631 {
632 GISMO_ENSURE( ( (i < 2) && (j < 2) && (k < 2) && (l < 2) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
633 GISMO_ENSURE(!comp,"Material model is not incompressible?");
634
635 return _Cijkl_impl<mat,imp>(i,j,k,l);
636 }
637
638 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
639 template <enum Material _mat, enum Implementation _imp>
640 typename std::enable_if<_mat==Material::SvK && _imp==Implementation::Analytical, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l) const641 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
642 {
643 // --------------------------
644 // Saint Venant Kirchhoff
645 // --------------------------
646 T lambda, mu, Cconstant;
647
648 mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
649 GISMO_ENSURE((1.-2.*m_parvals.at(1)) != 0, "Division by zero in construction of SvK material parameters! (1.-2.*m_parvals.at(1)) = "<<(1.-2.*m_parvals.at(1))<<"; m_parvals.at(1) = "<<m_parvals.at(1));
650 lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1))) ;
651 Cconstant = 2*lambda*mu/(lambda+2*mu);
652
653 return Cconstant*m_Acon_ori(i,j)*m_Acon_ori(k,l) + mu*(m_Acon_ori(i,k)*m_Acon_ori(j,l) + m_Acon_ori(i,l)*m_Acon_ori(j,k));
654 }
655
656 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
657 template <enum Material _mat, enum Implementation _imp>
658 typename std::enable_if<_mat==Material::NH && _imp==Implementation::Analytical, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l) const659 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
660 {
661 // --------------------------
662 // Neo-Hookean
663 // --------------------------
664 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
665 return mu*1./m_J0_sq*(2.*m_Gcon_def(i,j)*m_Gcon_def(k,l) + m_Gcon_def(i,k)*m_Gcon_def(j,l) + m_Gcon_def(i,l)*m_Gcon_def(j,k));
666 }
667
668 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
669 template <enum Material _mat, enum Implementation _imp>
670 typename std::enable_if<_mat==Material::MR && _imp==Implementation::Analytical, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l) const671 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
672 {
673 // --------------------------
674 // Mooney-Rivlin
675 // Parameter 3 is the ratio between c1 and c2.; c1 = m_parvals.at(2)*c2
676 // --------------------------
677 GISMO_ENSURE(m_numPars==3,"Mooney-Rivlin model needs to be a 3 parameter model");
678 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
679 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
680 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
681 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
682
683 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
684 T c2 = mu/(m_parvals.at(2) + 1);
685 T c1 = m_parvals.at(2)*c2;
686
687 T Gabcd = - 1./2. * ( m_Gcon_def(i,k)*m_Gcon_def(j,l) + m_Gcon_def(i,l)*m_Gcon_def(j,k) );
688
689 return (c1 + c2 * traceCt) *1./m_J0_sq*(2.*m_Gcon_def(i,j)*m_Gcon_def(k,l) + m_Gcon_def(i,k)*m_Gcon_def(j,l) + m_Gcon_def(i,l)*m_Gcon_def(j,k))// correct
690 - 2. * c2 / m_J0_sq * ( m_Gcon_ori(i,j) * m_Gcon_def(k,l) + m_Gcon_def(i,j)*m_Gcon_ori(k,l)) // correct
691 + 2. * c2 * m_J0_sq * ( Gabcd + m_Gcon_def(i,j)*m_Gcon_def(k,l) ); // Roohbakhshan
692 }
693
694 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
695 template <enum Material _mat, enum Implementation _imp>
696 typename std::enable_if<_imp==Implementation::Spectral, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l) const697 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
698 {
699 // --------------------------
700 // Stretch-based implementations
701 // --------------------------
702 T tmp = 0.0;
703 gsMatrix<T> C(3,3);
704 C.setZero();
705 C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
706 // C.block(0,0,2,2) = (m_gcov_def.transpose()*m_gcov_def).block(0,0,2,2);
707 // gsDebugVar(m_gcov_def.transpose()*m_gcov_def);
708 C(2,2) = 1./m_J0_sq;
709
710 this->_computeStretch(C);
711
712 tmp = 0.0;
713 for (index_t a = 0; a != 2; a++)
714 {
715 // C_iiii
716 tmp += _Cabcd(a,a,a,a)*(
717 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
718 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
719 );
720
721 for (index_t b = a+1; b != 2; b++)
722 {
723 // C_iijj = C_jjii
724 tmp += _Cabcd(a,a,b,b)*(
725 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
726 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
727 +
728 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
729 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
730 );
731
732 // C_ijij = Cjiji = Cijji = Cjiij
733 tmp += _Cabcd(a,b,a,b)*(
734 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
735 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
736 +
737 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
738 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
739 +
740 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
741 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
742 +
743 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
744 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
745 );
746 }
747 }
748 return tmp;
749 }
750
751 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
752 template <enum Material _mat, enum Implementation _imp>
753 typename std::enable_if<_imp==Implementation::Generalized, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l) const754 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
755 {
756 // --------------------------
757 // General implementations
758 // --------------------------
759 return 4.0 * _d2Psi(i,j,k,l) + 4.0 * _d2Psi(2,2,2,2)*math::pow(m_J0_sq,-2.0)*m_Gcon_def(i,j)*m_Gcon_def(k,l)
760 - 4.0/ m_J0_sq * ( _d2Psi(2,2,i,j)*m_Gcon_def(k,l) + _d2Psi(2,2,k,l)*m_Gcon_def(i,j) )
761 + 2.0 * _dPsi(2,2) / m_J0_sq * (2.*m_Gcon_def(i,j)*m_Gcon_def(k,l) + m_Gcon_def(i,k)*m_Gcon_def(j,l) + m_Gcon_def(i,l)*m_Gcon_def(j,k));
762 }
763
764 // template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
765 // template <enum Material _mat>
766 // // OTHER CASES!
767 // typename std::enable_if<(_mat >= 30) && !(_mat==0) && !(_mat==2) && !(_mat==3), T>::type
768 // gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
769 // {
770 // GISMO_ERROR("Material model unknown (model = "<<_mat<<"). Use gsMaterialMatrix<dim,T,matId,comp,mat,imp>::info() to see the options.");
771 // }
772
773 // else if (m_material==4)
774 // GISMO_ERROR("Material model 4 is not invariant-based! Use 14 instead...");
775 // else if (m_material==5)
776 // GISMO_ERROR("Material model 5 is only for compressible materials...");
777
778
779 // Condensation of the 3D tensor for compressible materials
780 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
781 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
782 {
783 GISMO_ENSURE(c.cols()==c.rows(),"Matrix c must be square");
784 GISMO_ENSURE(c.cols()==3,"Matrix c must be 3x3");
785 GISMO_ENSURE(cinv.cols()==cinv.rows(),"Matrix cinv must be square");
786 GISMO_ENSURE(cinv.cols()==3,"Matrix cinv must be 3x3");
787 GISMO_ENSURE( ( (i <2) && (j <2) && (k <2) && (l <2) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
788 GISMO_ENSURE(comp,"Material model is not compressible?");
789
790 return _Cijkl_impl<imp>(i,j,k,l,c,cinv);
791 }
792
793 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
794 template <enum Implementation _imp>
795 typename std::enable_if< _imp==Implementation::Spectral, T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const796 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
797 {
798 // static condensation is done before the projection
799 this->_computeStretch(c);
800 return _Cijkl3D(i,j,k,l,c,cinv);
801 }
802
803 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
804 template <enum Implementation _imp>
805 typename std::enable_if<!(_imp==Implementation::Spectral), T>::type
_Cijkl_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const806 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
807 {
808 return _Cijkl3D(i,j,k,l,c,cinv) - ( _Cijkl3D(i,j,2,2,c,cinv) * _Cijkl3D(2,2,k,l,c,cinv) ) / _Cijkl3D(2,2,2,2,c,cinv);
809 }
810
811 // 3D tensor for compressible materials
812 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
813 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
814 {
815 GISMO_ENSURE( ( (i <3) && (j <3) && (k <3) && (l <3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
816 return _Cijkl3D_impl<mat,imp>(i,j,k,l,c,cinv);
817 }
818
819 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
820 template <enum Material _mat, enum Implementation _imp>
821 typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const822 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
823 {
824 GISMO_ERROR("Compressible material matrix requested, but not needed. How?");
825 }
826
827 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
828 template <enum Material _mat, enum Implementation _imp>
829 typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const830 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
831 {
832 // --------------------------
833 // Neo-Hookean
834 // --------------------------
835
836 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
837 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
838
839 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
840 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
841 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
842 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
843 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
844 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
845 T I_1 = traceCt + c(2,2);
846 return 1.0 / 9.0 * mu * math::pow( m_J_sq , -1.0/3.0 ) * ( 2.0 * I_1 * ( cinv(i,j)*cinv(k,l) - 3.0 * dCinv )
847 - 6.0 *( m_Gcon_ori(i,j)*cinv(k,l) + cinv(i,j)*m_Gcon_ori(k,l) ) )
848 + K * ( m_J_sq*cinv(i,j)*cinv(k,l) + (m_J_sq-1)*dCinv );
849 }
850
851 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
852 template <enum Material _mat, enum Implementation _imp>
853 typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const854 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
855 {
856 // --------------------------
857 // Mooney-Rivlin
858 // --------------------------
859 GISMO_ENSURE(m_numPars==3,"Mooney-Rivlin model needs to be a 3 parameter model");
860
861 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
862 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
863
864 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
865 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
866 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
867 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
868 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
869 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
870 T I_1 = traceCt + c(2,2);
871 T I_2 = c(2,2) * traceCt + m_J0_sq;
872 T d2I_2 = idelta(i,2)*idelta(j,2)*idelta(k,2)*idelta(l,2)*( m_J0_sq*( cinv(i,j)*cinv(k,l) + dCinv ) )
873 + delta(i,2)*delta(j,2)*idelta(k,2)*idelta(l,2)*_dI_1(k,l)
874 + idelta(i,2)*idelta(j,2)*delta(k,2)*delta(l,2)*_dI_1(i,j);
875 T c2 = mu/(m_parvals.at(2) + 1);
876 T c1 = m_parvals.at(2)*c2;
877
878 return 1.0/9.0 * c1 * math::pow(m_J_sq, -1.0/3.0) * ( 2.0*I_1*cinv(i,j)*cinv(k,l) - 6.0*I_1*dCinv
879 - 6.0*_dI_1(i,j)*cinv(k,l) - 6.0*cinv(i,j)*_dI_1(k,l) ) // + 9*d2I_1 = 0
880 + 1.0/9.0 * c2 * math::pow(m_J_sq, -2.0/3.0) * ( 8.0*I_2*cinv(i,j)*cinv(k,l) - 12.0*I_2*dCinv
881 - 12.0*_dI_2(i,j,c,cinv)*cinv(k,l)- 12.0*cinv(i,j)*_dI_2(k,l,c,cinv)
882 + 18.0*d2I_2 )
883 + K * ( m_J_sq*cinv(i,j)*cinv(k,l) + (m_J_sq-1)*dCinv );
884 }
885
886 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
887 template <enum Material _mat, enum Implementation _imp>
888 typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const889 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
890 {
891 // --------------------------
892 // Neo-Hookean 2
893 // --------------------------
894 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
895 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
896 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
897
898 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
899 return - 2.0 * mu * dCinv + lambda * ( m_J_sq*cinv(i,j)*cinv(k,l) + (m_J_sq-1)*dCinv );
900 }
901
902 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
903 template <enum Material _mat, enum Implementation _imp>
904 typename std::enable_if<_imp == Implementation::Spectral, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const905 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
906 {
907 // --------------------------
908 // Stretch-based implementations
909 // --------------------------
910 T tmp;
911 if ( (i==2) && (j==2) && (k==2) && (l==2) ) // if C3333 (no static condensation)
912 tmp = _Cabcd(2,2,2,2);
913 else
914 {
915 tmp = 0.0;
916 T C = 0.0;
917 T C2222 = _Cabcd(2,2,2,2);
918 // T Cab22,C22ab;
919 for (index_t a = 0; a != 2; a++)
920 {
921 // C_iiii
922 // if (!((i==2 || j==2 || k==2 || l==2) && a!=2))
923 // C = _Cabcd(a,a,a,a);
924 C = _Cabcd(a,a,a,a) - math::pow(_Cabcd(2,2,a,a),2) / C2222;
925 tmp += C*(
926 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
927 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
928 );
929
930 for (index_t b = a+1; b != 2; b++)
931 {
932 // C_iijj
933 // C = _Cabcd(a,a,b,b);
934 C = _Cabcd(a,a,b,b) - _Cabcd(a,a,2,2) * _Cabcd(2,2,b,b) / C2222;
935 tmp += C*(
936 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
937 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
938 +
939 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
940 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
941 );
942
943 // C_ijij = Cjiji = Cijji = Cjiij
944 // C = _Cabcd(a,b,a,b);
945 C = _Cabcd(a,b,a,b) - math::pow(_Cabcd(2,2,a,b),2) / C2222;
946 tmp += C*(
947 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
948 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
949 +
950 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
951 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
952 +
953 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(b)) )*
954 ( m_gcon_ori.col(k).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(a)) )
955 +
956 ( m_gcon_ori.col(i).dot(m_stretchvec.col(b)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )*
957 ( m_gcon_ori.col(k).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(l).dot(m_stretchvec.col(b)) )
958 );
959 }
960 }
961 }
962
963 return tmp;
964 }
965
966 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
967 template <enum Material _mat, enum Implementation _imp>
968 typename std::enable_if<_imp == Implementation::Generalized, T>::type
_Cijkl3D_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const969 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cijkl3D_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
970 {
971 // --------------------------
972 // General implementations
973 // --------------------------
974 return 4.0 * _d2Psi(i,j,k,l,c,cinv);
975 }
976
977 /*
978 Available class members:
979 - m_parvals.at(0)
980 - m_parvals.at(1)
981 - m_metric
982 - m_metric_def
983 - m_J0
984 - m_J
985 - m_Cinv
986 */
987 // template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
988 // T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j) const { _Sij(i,j,NULL,NULL); }
989
990 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
991 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j) const
992 {
993 return _Sij_impl<mat,imp>(i,j);
994 }
995
996 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
997 template <enum Material _mat, enum Implementation _imp>
998 typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j) const999 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1000 {
1001 GISMO_ERROR("Not implemented");
1002 }
1003
1004 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1005 template <enum Material _mat, enum Implementation _imp>
1006 typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j) const1007 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1008 {
1009 // --------------------------
1010 // Neo-Hoookean
1011 // --------------------------
1012 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
1013 return mu * (m_Gcon_ori(i,j) - 1./m_J0_sq * m_Gcon_def(i,j) );
1014 }
1015
1016 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1017 template <enum Material _mat, enum Implementation _imp>
1018 typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j) const1019 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1020 {
1021 // --------------------------
1022 // Mooney-Rivlin
1023 // Parameter 3 is the ratio between c1 and c2.
1024 // --------------------------
1025 GISMO_ENSURE(m_numPars==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1026 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
1027 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1028 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1029 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1030 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1031
1032 T c2 = mu/(m_parvals.at(2)+1);
1033 T c1 = m_parvals.at(2)*c2;
1034
1035 return c1 * ( m_Gcon_ori(i,j) - 1/m_J0_sq * m_Gcon_def(i,j) )
1036 + c2 / m_J0_sq * (m_Gcon_ori(i,j) - traceCt * m_Gcon_def(i,j) ) + c2 * m_J0_sq * m_Gcon_def(i,j);
1037 }
1038
1039 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1040 template <enum Material _mat, enum Implementation _imp>
1041 typename std::enable_if<_imp == Implementation::Spectral, T>::type
_Sij_impl(const index_t i,const index_t j) const1042 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1043 {
1044 T tmp = 0.0;
1045 gsMatrix<T> C(3,3);
1046 C.setZero();
1047 C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
1048 C(2,2) = 1./m_J0_sq;
1049
1050 this->_computeStretch(C);
1051
1052 for (index_t a = 0; a != 2; a++)
1053 {
1054 tmp += _Sa(a)*(
1055 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )
1056 );
1057 }
1058 return tmp;
1059 }
1060
1061 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1062 template <enum Material _mat, enum Implementation _imp>
1063 typename std::enable_if<_imp == Implementation::Generalized, T>::type
_Sij_impl(const index_t i,const index_t j) const1064 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1065 {
1066 // --------------------------
1067 // Generalized
1068 // --------------------------
1069 return 2.0 * _dPsi(i,j) - 2.0 * _dPsi(2,2) * math::pow(m_J0_sq,-1.0)*m_Gcon_def(i,j);
1070 }
1071
1072 //--------------------------------------------------------------------------------------------------------------------------------------
1073
1074 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1075 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1076 {
1077 return _Sij_impl<mat,imp>(i,j,c,cinv);
1078 }
1079
1080 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1081 template <enum Material _mat, enum Implementation _imp>
1082 typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1083 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1084 {
1085 // --------------------------
1086 // Neo-Hoookean
1087 // --------------------------
1088 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
1089 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1090 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1091 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1092 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1093 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1094 T I_1 = traceCt + c(2,2);
1095
1096 return mu * math::pow( m_J_sq , -1.0/3.0 ) * ( m_Gcon_ori(i,j) - 1.0/3.0 * I_1 * cinv(i,j) )
1097 + K * 0.5 * ( m_J_sq - 1.0 ) * cinv(i,j);
1098 }
1099
1100 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1101 template <enum Material _mat, enum Implementation _imp>
1102 typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1103 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1104 {
1105 // --------------------------
1106 // Mooney-Rivlin
1107 // Parameter 3 is the ratio between c1 and c2.
1108 // --------------------------
1109 GISMO_ENSURE(m_numPars==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1110 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
1111 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1112 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1113 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1114 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1115 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1116 T I_1 = traceCt + c(2,2);
1117 T I_2 = c(2,2) * traceCt + m_J0_sq;
1118
1119 T c2 = mu/(m_parvals.at(2)+1);
1120 T c1 = m_parvals.at(2)*c2;
1121
1122 return c1 * math::pow( m_J_sq , -1.0/3.0 ) * ( m_Gcon_ori(i,j) - 1.0/3.0 * I_1 * cinv(i,j) )
1123 + c2 * math::pow( m_J_sq , -2.0/3.0 ) * ( _dI_2(i,j,c,cinv)- 2.0/3.0 * I_2 * cinv(i,j) )
1124 + K * 0.5 * ( m_J_sq - 1.0 ) * cinv(i,j);
1125 }
1126
1127 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1128 template <enum Material _mat, enum Implementation _imp>
1129 typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
_Sij_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1130 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1131 {
1132 // --------------------------
1133 // Neo-Hookean 2
1134 // --------------------------
1135 T mu = m_parvals.at(0) / (2.*(1. + m_parvals.at(1)));
1136 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
1137 return mu * m_Gcon_ori(i,j) - mu * cinv(i,j) + lambda / 2.0 * ( m_J_sq - 1 ) * cinv(i,j);
1138 }
1139
1140 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1141 template <enum Material _mat, enum Implementation _imp>
1142 typename std::enable_if<_imp == Implementation::Spectral, T>::type
_Sij_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1143 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1144 {
1145 T tmp = 0.0;
1146 this->_computeStretch(c);
1147 for (index_t a = 0; a != 3; a++)
1148 {
1149 tmp += _Sa(a)*(
1150 ( m_gcon_ori.col(i).dot(m_stretchvec.col(a)) )*( m_gcon_ori.col(j).dot(m_stretchvec.col(a)) )
1151 );
1152 }
1153 return tmp;
1154 }
1155
1156 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1157 template <enum Material _mat, enum Implementation _imp>
1158 typename std::enable_if<_imp == Implementation::Generalized, T>::type
_Sij_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1159 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1160 {
1161 // --------------------------
1162 // Generalized
1163 // --------------------------
1164 return 2.0 * _dPsi(i,j,c,cinv);
1165 }
1166
1167 //--------------------------------------------------------------------------------------------------------------------------------------
1168
1169 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1170 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sii(const index_t i) const // principle stresses
1171 {
1172 gsMatrix<T> C(3,3);
1173 C.setZero();
1174 C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
1175 C(2,2) = 1./m_J0_sq;
1176 return _Sii(i,C);
1177 }
1178
1179 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1180 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sii(const index_t i, const gsMatrix<T> & c) const
1181 {
1182 this->_computeStretch(c);
1183 return _Sa(i);
1184 }
1185
1186 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Compressible_matrix(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const1187 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Compressible_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1188 {
1189 // Input: j index in-plane point
1190 // z out-of-plane coordinate (through thickness) in R1 (z)
1191 // Output: (n=u.cols(), m=z.rows())
1192 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1193 this->_computePoints(patch,u);
1194 gsMatrix<T> result(9, u.cols() * z.rows());
1195 result.setZero();
1196
1197 for (index_t k=0; k!=u.cols(); k++)
1198 {
1199 // Evaluate material properties on the quadrature point
1200 for (index_t v=0; v!=m_parmat.rows(); v++)
1201 m_parvals.at(v) = m_parmat(v,k);
1202
1203 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1204 {
1205 // this->computeMetric(i,z.at(j),true,true); // on point i, on height z(0,j)
1206 this->_getMetric(k, z(j, k) * m_Tmat(0, k)); // on point i, on height z(0,j)
1207
1208 // Define objects
1209 gsMatrix<T,3,3> c, cinv;
1210 T S33, C3333, dc33;
1211 // T S33_old;
1212 S33 = 0.0;
1213 dc33 = 0.0;
1214 C3333 = 1.0;
1215
1216 index_t itmax = 100;
1217 T tol = 1e-10;
1218
1219 // Initialize c
1220 c.setZero();
1221 c.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
1222 c(2,2) = math::pow(m_J0_sq,-1.0); // c33
1223 // c(2,2) = 1.0; // c33
1224 cinv.setZero();
1225 cinv.block(0,0,2,2) = m_Gcon_def.block(0,0,2,2);
1226 cinv(2,2) = 1.0/c(2,2);
1227
1228 m_J_sq = m_J0_sq * c(2,2);
1229 S33 = _Sij(2,2,c,cinv);
1230 // S33_old = (S33 == 0.0) ? 1.0 : S33;
1231 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
1232
1233 dc33 = -2. * S33 / C3333;
1234 for (index_t it = 0; it < itmax; it++)
1235 {
1236 c(2,2) += dc33;
1237
1238 GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
1239 cinv(2,2) = 1.0/c(2,2);
1240
1241 m_J_sq = m_J0_sq * c(2,2) ;
1242
1243 S33 = _Sij(2,2,c,cinv);
1244 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
1245
1246 dc33 = -2. * S33 / C3333;
1247 // if (abs(S33/S33_old) < tol)
1248 if (abs(dc33) < tol)
1249 {
1250 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(j*u.cols()+k,3,3);
1251 /*
1252 C = C1111, C1122, C1112
1253 symm, C2222, C2212
1254 symm, symm, C1212
1255 */
1256 C(0,0) = _Cijkl(0,0,0,0,c,cinv); // C1111
1257 C(1,1) = _Cijkl(1,1,1,1,c,cinv); // C2222
1258 C(2,2) = _Cijkl(0,1,0,1,c,cinv); // C1212
1259 C(1,0) = C(0,1) = _Cijkl(0,0,1,1,c,cinv); // C1122
1260 C(2,0) = C(0,2) = _Cijkl(0,0,0,1,c,cinv); // C1112
1261 C(2,1) = C(1,2) = _Cijkl(1,1,0,1,c,cinv); // C2212
1262
1263 break;
1264 }
1265 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<" and tolerance = "<<tol<<"\n");
1266 }
1267 }
1268 }
1269 return result;
1270 }
1271
1272 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Compressible_vector(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const1273 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Compressible_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1274 {
1275 // Input: j index in-plane point
1276 // z out-of-plane coordinate (through thickness) in R1 (z)
1277 // Output: (n=u.cols(), m=z.rows())
1278 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1279 this->_computePoints(patch,u);
1280 gsMatrix<T> result(3, u.cols() * z.rows());
1281 result.setZero();
1282
1283 for (index_t k=0; k!=u.cols(); k++)
1284 {
1285 // Evaluate material properties on the quadrature point
1286 for (index_t v=0; v!=m_parmat.rows(); v++)
1287 m_parvals.at(v) = m_parmat(v,k);
1288
1289 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1290 {
1291 // this->computeMetric(i,z.at(j),true,true); // on point i, on height z(0,j)
1292 this->_getMetric(k, z(j, k) * m_Tmat(0, k)); // on point i, on height z(0,j)
1293
1294 // Define objects
1295 gsMatrix<T,3,3> c, cinv;
1296 T S33, C3333, dc33;
1297 // T S33_old;
1298 S33 = 0.0;
1299 dc33 = 0.0;
1300 C3333 = 1.0;
1301
1302 index_t itmax = 100;
1303 T tol = 1e-10;
1304
1305 // Initialize c
1306 c.setZero();
1307 c.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
1308 c(2,2) = math::pow(m_J0_sq,-1.0); // c33
1309 // c(2,2) = 1.0; // c33
1310 cinv.setZero();
1311 cinv.block(0,0,2,2) = m_Gcon_def.block(0,0,2,2);
1312 cinv(2,2) = 1.0/c(2,2);
1313
1314 m_J_sq = m_J0_sq * c(2,2);
1315 S33 = _Sij(2,2,c,cinv);
1316 // S33_old = (S33 == 0.0) ? 1.0 : S33;
1317 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
1318
1319 dc33 = -2. * S33 / C3333;
1320 for (index_t it = 0; it < itmax; it++)
1321 {
1322 c(2,2) += dc33;
1323
1324 GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
1325 cinv(2,2) = 1.0/c(2,2);
1326
1327 m_J_sq = m_J0_sq * c(2,2) ;
1328
1329 S33 = _Sij(2,2,c,cinv);
1330 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
1331
1332 dc33 = -2. * S33 / C3333;
1333 // if (abs(S33/S33_old) < tol)
1334 if (abs(dc33) < tol)
1335 {
1336
1337 result(0,j*u.cols()+k) = _Sij(0,0,c,cinv); // S11
1338 result(1,j*u.cols()+k) = _Sij(1,1,c,cinv); // S22
1339 result(2,j*u.cols()+k) = _Sij(0,1,c,cinv); // S12
1340
1341 break;
1342 }
1343 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<" and tolerance = "<<tol<<"\n");
1344 }
1345 }
1346 }
1347 return result;
1348 }
1349
1350 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
_eval_Compressible_pstress(const index_t patch,const gsMatrix<T> & u,const gsMatrix<T> & z) const1351 gsMatrix<T> gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_eval_Compressible_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1352 {
1353 // Input: j index in-plane point
1354 // z out-of-plane coordinate (through thickness) in R1 (z)
1355 // Output: (n=u.cols(), m=z.rows())
1356 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1357 this->_computePoints(patch,u);
1358 gsMatrix<T> result(2, u.cols() * z.rows());
1359 result.setZero();
1360
1361 for (index_t k=0; k!=u.cols(); k++)
1362 {
1363 // Evaluate material properties on the quadrature point
1364 for (index_t v=0; v!=m_parmat.rows(); v++)
1365 m_parvals.at(v) = m_parmat(v,k);
1366
1367 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1368 {
1369 // this->computeMetric(i,z.at(j),true,true); // on point i, on height z(0,j)
1370 this->_getMetric(k, z(j, k) * m_Tmat(0, k)); // on point i, on height z(0,j)
1371
1372 // Define objects
1373 gsMatrix<T,3,3> c, cinv;
1374 T S33, C3333, dc33;
1375 // T S33_old;
1376 S33 = 0.0;
1377 dc33 = 0.0;
1378 C3333 = 1.0;
1379
1380 index_t itmax = 100;
1381 T tol = 1e-10;
1382
1383 // Initialize c
1384 c.setZero();
1385 c.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
1386 c(2,2) = math::pow(m_J0_sq,-1.0); // c33
1387 // c(2,2) = 1.0; // c33
1388
1389 cinv.setZero();
1390 cinv.block(0,0,2,2) = m_Gcon_def.block(0,0,2,2);
1391 cinv(2,2) = 1.0/c(2,2);
1392
1393 m_J_sq = m_J0_sq * c(2,2);
1394 S33 = _Sij(2,2,c,cinv);
1395 // S33_old = (S33 == 0.0) ? 1.0 : S33;
1396 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
1397
1398 dc33 = -2. * S33 / C3333;
1399 for (index_t it = 0; it < itmax; it++)
1400 {
1401 c(2,2) += dc33;
1402
1403 GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
1404 cinv(2,2) = 1.0/c(2,2);
1405
1406 m_J_sq = m_J0_sq * c(2,2) ;
1407
1408 S33 = _Sij(2,2,c,cinv);
1409 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
1410
1411 dc33 = -2. * S33 / C3333;
1412 // if (abs(S33/S33_old) < tol)
1413 if (abs(dc33) < tol)
1414 {
1415 GISMO_ENSURE(imp==Implementation::Spectral, "Only available for stretch-based materials.");
1416
1417 result(0,j*u.cols()+k) = _Sii(0,c); // S11
1418 result(1,j*u.cols()+k) = _Sii(1,c); // S22
1419
1420 break;
1421 }
1422 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<" and tolerance = "<<tol<<"\n");
1423 }
1424 }
1425 }
1426 return result;
1427 }
1428
1429 // ---------------------------------------------------------------------------------------------------------------------------------
1430 // INCOMPRESSIBLE
1431 // ---------------------------------------------------------------------------------------------------------------------------------
1432 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1433 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi(const index_t i, const index_t j) const
1434 {
1435 GISMO_ENSURE( ( (i < 3) && (j < 3) ) , "Index out of range. i="<<i<<", j="<<j);
1436 GISMO_ENSURE(!comp,"Material model is not incompressible?");
1437 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
1438 return _dPsi_impl<mat>(i,j);
1439 }
1440
1441 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1442 template <enum Material _mat>
1443 typename std::enable_if<_mat==Material::NH, T>::type
_dPsi_impl(const index_t i,const index_t j) const1444 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j) const
1445 {
1446 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1447 return 0.5 * mu * m_Gcon_ori(i,j);
1448 }
1449
1450 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1451 template <enum Material _mat>
1452 typename std::enable_if<_mat==Material::MR, T>::type
_dPsi_impl(const index_t i,const index_t j) const1453 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j) const
1454 {
1455 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1456 T c2 = mu/(m_parvals.at(2)+1);
1457 T c1 = m_parvals.at(2)*c2;
1458 T tmp;
1459 if ((i==2) && (j==2))
1460 {
1461 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1462 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1463 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1464 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1465 tmp = c1/2.0 + c2 / 2.0 * traceCt;
1466 }
1467 else
1468 tmp = c1 / 2. * m_Gcon_ori(i,j) + c2 / 2. * ( 1. / m_J0_sq * m_Gcon_ori(i,j) + m_J0_sq * m_Gcon_def(i,j) );
1469
1470 return tmp;
1471 }
1472
1473
1474 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1475 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi(const index_t i, const index_t j, const index_t k, const index_t l) const
1476 {
1477 GISMO_ENSURE( ( (i < 3) && (j < 3) && (k < 3) && (l < 3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
1478 GISMO_ENSURE(!comp,"Material model is not incompressible?");
1479 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
1480 return _d2Psi_impl<mat>(i,j,k,l);
1481 }
1482
1483 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1484 template <enum Material _mat>
1485 typename std::enable_if<_mat==Material::NH, T>::type
_d2Psi_impl(const index_t i,const index_t j,const index_t k,const index_t l) const1486 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
1487 {
1488 return 0.0;
1489 }
1490
1491 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1492 template <enum Material _mat>
1493 typename std::enable_if<_mat==Material::MR, T>::type
_d2Psi_impl(const index_t i,const index_t j,const index_t k,const index_t l) const1494 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_impl(const index_t i, const index_t j, const index_t k, const index_t l) const
1495 {
1496 T tmp;
1497 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1498 T c2 = mu/(m_parvals.at(2)+1);
1499 // T c1 = m_parvals.at(2)*c2;
1500 if ( ((i==2) && (j==2)) && !((k==2) || (l==2)) ) // _dPsi/d22dkl
1501 tmp = c2 / 2.0 * m_Gcon_ori(k,l);
1502 else if ( !((i==2) && (j==2)) && ((k==2) || (l==2)) ) // _dPsi/dijd22
1503 tmp = c2 / 2.0 * m_Gcon_ori(i,j);
1504 else if ( ((i==2) && (j==2)) && ((k==2) || (l==2)) ) // _dPsi/d22d22
1505 tmp = 0.0;
1506 else
1507 {
1508 T Gabcd = - 1./2. * ( m_Gcon_def(i,k)*m_Gcon_def(j,l) + m_Gcon_def(i,l)*m_Gcon_def(j,k) );
1509 tmp = c2 / 2.0 * m_J0_sq * ( Gabcd + m_Gcon_def(i,j)*m_Gcon_def(k,l) );
1510 }
1511 return tmp;
1512 }
1513
1514 // ---------------------------------------------------------------------------------------------------------------------------------
1515 // COMPRESSIBLE
1516 // ---------------------------------------------------------------------------------------------------------------------------------
1517
1518 //--------------------------------------------------------------------------------------------------------------------------------------
1519
1520 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1521 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dI_1(const index_t i, const index_t j) const
1522 {
1523 return m_Gcon_ori(i,j);
1524 }
1525
1526 //--------------------------------------------------------------------------------------------------------------------------------------
1527
1528 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1529 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dI_2(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1530 {
1531 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1532 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1533 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1534 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1535 return idelta(i,2)*idelta(j,2)*( c(2,2)*_dI_1(i,j) + m_J0_sq*cinv(i,j) ) + delta(i,2)*delta(j,2)*traceCt;
1536 }
1537
1538 //--------------------------------------------------------------------------------------------------------------------------------------
1539
1540 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1541 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1542 {
1543 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
1544 return _dPsi_impl<mat>(i,j,c,cinv);
1545 }
1546
1547 //--------------------------------------------------------------------------------------------------------------------------------------
1548
1549 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1550 template <enum Material _mat>
1551 typename std::enable_if<_mat==Material::NH, T>::type
_dPsi_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1552 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1553 {
1554 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1555 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1556 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1557 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1558 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1559 T I_1 = traceCt + c(2,2);
1560 return mu/2.0 * math::pow(m_J_sq,-1./3.) * ( - 1.0/3.0 * I_1 * cinv(i,j) + _dI_1(i,j) ) + _dPsi_vol(i,j,c,cinv);
1561 }
1562
1563 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1564 template <enum Material _mat>
1565 typename std::enable_if<_mat==Material::MR, T>::type
_dPsi_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1566 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1567 {
1568 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1569 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1570 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1571 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1572 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1573 T I_1 = traceCt + c(2,2);
1574 T I_2 = c(2,2) * traceCt + m_J0_sq;
1575 T c2= mu/(m_parvals.at(2)+1);
1576 T c1= m_parvals.at(2)*c2;
1577 return c1/2.0 * math::pow(m_J_sq,-1./3.) * ( - 1.0/3.0 * I_1 * cinv(i,j) + _dI_1(i,j) )
1578 + c2/2.0 * math::pow(m_J_sq,-2./3.) * ( - 2.0/3.0 * I_2 * cinv(i,j) + _dI_2(i,j,c,cinv) )
1579 + _dPsi_vol(i,j,c,cinv);
1580 }
1581
1582 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1583 template <enum Material _mat>
1584 typename std::enable_if<_mat==Material::NH_ext, T>::type
_dPsi_impl(const index_t i,const index_t j,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1585 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1586 {
1587 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1588 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
1589 return mu / 2.0 * _dI_1(i,j) - mu / 2.0 * cinv(i,j) + lambda / 4.0 * ( m_J_sq - 1 ) * cinv(i,j);
1590 }
1591
1592 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1593 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_vol(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1594 {
1595 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1596 return K * 0.25 * (m_J_sq - 1.0) * cinv(i,j);
1597 }
1598
1599 //--------------------------------------------------------------------------------------------------------------------------------------
1600
1601 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1602 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1603 {
1604 GISMO_ENSURE( ( (i < 3) && (j < 3) && (k < 3) && (l < 3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
1605 GISMO_ENSURE(comp,"Material model is not compressible?");
1606 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
1607 return _d2Psi_impl<mat>(i,j,k,l,c,cinv);
1608 }
1609
1610 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1611 template <enum Material _mat>
1612 typename std::enable_if<_mat==Material::NH, T>::type
_d2Psi_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1613 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1614 {
1615 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1616 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1617 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1618 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1619 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1620 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1621 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1622 T I_1 = traceCt + c(2,2);
1623 return 1.0/9.0 * mu / 2.0 * math::pow(m_J_sq, -1.0/3.0) * ( I_1*cinv(i,j)*cinv(k,l)
1624 - 3.0*_dI_1(i,j)*cinv(k,l) - 3.0*cinv(i,j)*_dI_1(k,l)
1625 - 3.0*I_1*dCinv )
1626 + _d2Psi_vol(i,j,k,l,c,cinv);
1627 }
1628
1629 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1630 template <enum Material _mat>
1631 typename std::enable_if<_mat==Material::MR, T>::type
_d2Psi_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1632 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1633 {
1634 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1635 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1636 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1637 T traceCt = m_Gcov_def(0,0)*m_Gcon_ori(0,0) +
1638 m_Gcov_def(0,1)*m_Gcon_ori(0,1) +
1639 m_Gcov_def(1,0)*m_Gcon_ori(1,0) +
1640 m_Gcov_def(1,1)*m_Gcon_ori(1,1);
1641 T I_1 = traceCt + c(2,2);
1642 T I_2 = c(2,2) * traceCt + m_J0_sq;
1643 T d2I_2 = idelta(i,2)*idelta(j,2)*idelta(k,2)*idelta(l,2)*( m_J0_sq*( cinv(i,j)*cinv(k,l) + dCinv ) )
1644 + delta(i,2)*delta(j,2)*idelta(k,2)*idelta(l,2)*_dI_1(k,l)
1645 + idelta(i,2)*idelta(j,2)*delta(k,2)*delta(l,2)*_dI_1(i,j);
1646 T c2 = mu/(m_parvals.at(2)+1);
1647 T c1 = m_parvals.at(2)*c2;
1648 // c1 = 0;
1649 return
1650 1.0/9.0 * c1 / 2.0 * math::pow(m_J_sq, -1.0/3.0) * ( I_1*cinv(i,j)*cinv(k,l)
1651 - 3.0*_dI_1(i,j)*cinv(k,l) - 3.0*cinv(i,j)*_dI_1(k,l)
1652 - 3.0*I_1*dCinv ) // + 9*d2I_1 = 0
1653 + 1.0/9.0 * c2 / 2.0 * math::pow(m_J_sq, -2.0/3.0) * ( 4.0*I_2*cinv(i,j)*cinv(k,l) - 6.0*I_2*dCinv
1654 - 6.0*_dI_2(i,j,c,cinv)*cinv(k,l)- 6.0*cinv(i,j)*_dI_2(k,l,c,cinv)
1655 + 9.0*d2I_2 )
1656 + _d2Psi_vol(i,j,k,l,c,cinv);
1657 }
1658
1659 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1660 template <enum Material _mat>
1661 typename std::enable_if<_mat==Material::NH_ext, T>::type
_d2Psi_impl(const index_t i,const index_t j,const index_t k,const index_t l,const gsMatrix<T> & c,const gsMatrix<T> & cinv) const1662 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_impl(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1663 {
1664 T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
1665 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1666 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
1667 return - mu / 2.0 * dCinv + lambda / 4.0 * ( m_J_sq*cinv(i,j)*cinv(k,l) + (m_J_sq-1.0)*dCinv );
1668 }
1669
1670 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1671 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_vol(const index_t i, const index_t j, const index_t k, const index_t l, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1672 {
1673 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1674 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1675 return K * 0.25 * ( m_J_sq*cinv(i,j)*cinv(k,l) + (m_J_sq-1.0)*dCinv );
1676 }
1677
1678 // ---------------------------------------------------------------------------------------------------------------------------------
1679 // STRETCHES
1680 // ---------------------------------------------------------------------------------------------------------------------------------
1681
1682 //--------------------------------------------------------------------------------------------------------------------------------------
1683
1684 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1685 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da(const index_t a) const
1686 {
1687 GISMO_ENSURE( a < 3 , "Index out of range. a="<<a);
1688 return _dPsi_da_impl<mat,comp>(a);
1689 }
1690
1691 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1692 template <enum Material _mat, bool _comp>
1693 typename std::enable_if<_comp && (_mat==Material::NH), T>::type
_dPsi_da_impl(const index_t a) const1694 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1695 {
1696 GISMO_ENSURE( a < 3 , "Index out of range. a="<<a);
1697 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1698 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1699 T _dI_1a = 2*m_stretches(a);
1700
1701 return mu/2.0 * math::pow(m_J_sq,-1./3.) * ( -2./3. * I_1 / m_stretches(a) + _dI_1a )
1702 + _dPsi_da_vol(a);
1703 }
1704
1705 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1706 template <enum Material _mat, bool _comp>
1707 typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
_dPsi_da_impl(const index_t a) const1708 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1709 {
1710 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1711 T _dI_1a = 2*m_stretches(a);
1712 return mu/2 * _dI_1a;
1713 }
1714
1715 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1716 template <enum Material _mat, bool _comp>
1717 typename std::enable_if<_comp && (_mat==Material::MR), T>::type
_dPsi_da_impl(const index_t a) const1718 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1719 {
1720 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1721 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1722 T _dI_1a = 2*m_stretches(a);
1723 T I_2 = math::pow(m_stretches(0),2.)*math::pow(m_stretches(1),2.)
1724 + math::pow(m_stretches(1),2.)*math::pow(m_stretches(2),2.)
1725 + math::pow(m_stretches(0),2.)*math::pow(m_stretches(2),2.);
1726 T _dI_2a = 2*m_stretches(a)*( I_1 - math::pow(m_stretches(a),2.0) );
1727
1728 T c2= mu/(m_parvals.at(2)+1);
1729 T c1= m_parvals.at(2)*c2;
1730 return c1/2.0 * math::pow(m_J_sq,-1./3.) * ( -2./3. * I_1 / m_stretches(a) + _dI_1a )
1731 + c2/2.0 * math::pow(m_J_sq,-2./3.) * ( -4./3. * I_2 / m_stretches(a) + _dI_2a )
1732 + _dPsi_da_vol(a);
1733 }
1734
1735 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1736 template <enum Material _mat, bool _comp>
1737 typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
_dPsi_da_impl(const index_t a) const1738 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1739 {
1740 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1741 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1742 T _dI_1a = 2*m_stretches(a);
1743 T _dI_2a = 2*m_stretches(a)*( I_1 - math::pow(m_stretches(a),2.0) );
1744
1745 T c2 = mu/(m_parvals.at(2)+1);
1746 T c1 = m_parvals.at(2)*c2;
1747 return c1/2.0*_dI_1a + c2/2.0*_dI_2a;
1748 }
1749
1750 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1751 template <enum Material _mat, bool _comp>
1752 typename std::enable_if<_comp && (_mat==Material::OG), T>::type
_dPsi_da_impl(const index_t a) const1753 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1754 {
1755 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1756 T tmp = 0.0;
1757 int n = (m_numPars-2)/2;
1758 T alpha_i, mu_i, Lambda;
1759 for (index_t k=0; k!=n; k++)
1760 {
1761 alpha_i = m_parvals.at(2*(k+1)+1);
1762 mu_i = m_parvals.at(2*(k+1));
1763 Lambda = math::pow(m_stretches(0),alpha_i) + math::pow(m_stretches(1),alpha_i) + math::pow(m_stretches(2),alpha_i);
1764 tmp += mu_i * math::pow(m_J_sq,-alpha_i/6.0) * ( math::pow(m_stretches(a),alpha_i-1) - 1./3. * 1./m_stretches(a) * Lambda );
1765 }
1766 return tmp + _dPsi_da_vol(a);
1767 }
1768 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1769 template <enum Material _mat, bool _comp>
1770 typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
_dPsi_da_impl(const index_t a) const1771 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1772 {
1773 T tmp = 0.0;
1774 int n = (m_numPars-2)/2;
1775 T alpha_i, mu_i;
1776 for (index_t k=0; k!=n; k++)
1777 {
1778 alpha_i = m_parvals.at(2*(k+1)+1);
1779 mu_i = m_parvals.at(2*(k+1));
1780 tmp += mu_i*math::pow(m_stretches(a),alpha_i-1);
1781 }
1782 return tmp;
1783 }
1784
1785 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1786 template <enum Material _mat, bool _comp>
1787 typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
_dPsi_da_impl(const index_t a) const1788 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
1789 {
1790 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1791 T _dI_1a = 2*m_stretches(a);
1792 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1793 // choose compressibility function (and parameter)
1794 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
1795
1796 return mu/2.0 * _dI_1a - mu / m_stretches(a) + lambda / (m_stretches(a)*2) * (m_J_sq-1.0);
1797 }
1798
1799 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1800 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi_da_vol(const index_t a) const
1801 {
1802 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1803 T beta = -2.0;
1804 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1805 return K / (m_stretches(a)*beta) * (1.0 - math::pow(m_J_sq,-beta/2.0));
1806 }
1807
1808 //--------------------------------------------------------------------------------------------------------------------------------------
1809
1810 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1811 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab(const index_t a, const index_t b) const
1812 {
1813 GISMO_ENSURE( ( (a < 3) && (b < 3) ) , "Index out of range. a="<<a<<", b="<<b);
1814 return _d2Psi_dab_impl<mat,comp>(a,b);
1815 }
1816
1817 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1818 template <enum Material _mat, bool _comp>
1819 typename std::enable_if<_comp && (_mat==Material::NH), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1820 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1821 {
1822 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1823
1824 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1825 T _dI_1a = 2*m_stretches(a);
1826 T _dI_1b = 2*m_stretches(b);
1827 T d2I_1 = 2*delta(a,b);
1828 return mu/2.0 * math::pow(m_J_sq,-1./3.) * (
1829 -2./3. * 1. / m_stretches(b) * ( -2./3. * I_1 / m_stretches(a) + _dI_1a )
1830 -2./3. * 1. / m_stretches(a) * _dI_1b
1831 +d2I_1
1832 +2./3. * delta(a,b) * I_1 / (m_stretches(a)*m_stretches(a))
1833 )
1834 + _d2Psi_dab_vol(a,b);
1835 }
1836
1837 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1838 template <enum Material _mat, bool _comp>
1839 typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1840 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1841 {
1842 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1843 T d2I_1 = 2*delta(a,b);
1844 return mu/2 * d2I_1;
1845 }
1846
1847 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1848 template <enum Material _mat, bool _comp>
1849 typename std::enable_if<_comp && (_mat==Material::MR), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1850 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1851 {
1852 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1853
1854 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1855 T _dI_1a = 2*m_stretches(a);
1856 T _dI_1b = 2*m_stretches(b);
1857 T d2I_1 = 2*delta(a,b);
1858
1859 T I_2 = math::pow(m_stretches(0),2.)*math::pow(m_stretches(1),2.)
1860 + math::pow(m_stretches(1),2.)*math::pow(m_stretches(2),2.)
1861 + math::pow(m_stretches(0),2.)*math::pow(m_stretches(2),2.);
1862 T _dI_2a = 2*m_stretches(a)*( I_1 - math::pow(m_stretches(a),2.0) );
1863 T _dI_2b = 2*m_stretches(b)*( I_1 - math::pow(m_stretches(b),2.0) );
1864 T d2I_2 = idelta(a,b)*4.0*m_stretches(a)*m_stretches(b) + delta(a,b)*2.0*(I_1 - m_stretches(a)*m_stretches(a));
1865
1866 T c2 = mu/(m_parvals.at(2)+1);
1867 T c1 = m_parvals.at(2)*c2;
1868 return
1869 c1/2.0 * math::pow(m_J_sq,-1./3.) * (
1870 -2./3. * 1. / m_stretches(b) * ( -2./3. * I_1 / m_stretches(a) + _dI_1a )
1871 -2./3. * 1. / m_stretches(a) * _dI_1b
1872 +d2I_1
1873 +2./3. * delta(a,b) * I_1 / (m_stretches(a)*m_stretches(a))
1874 )
1875 + c2/2.0 * math::pow(m_J_sq,-2./3.) * (
1876 -4./3. * 1. / m_stretches(b) * ( -4./3. * I_2 / m_stretches(a) + _dI_2a )
1877 -4./3. * 1. / m_stretches(a) * _dI_2b
1878 +d2I_2
1879 +4./3. * delta(a,b) * I_2 / (m_stretches(a)*m_stretches(a))
1880 )
1881 + _d2Psi_dab_vol(a,b);
1882 }
1883
1884 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1885 template <enum Material _mat, bool _comp>
1886 typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1887 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1888 {
1889 GISMO_ENSURE(m_numPars==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1890 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1891
1892 T I_1 = m_stretches(0)*m_stretches(0) + m_stretches(1)*m_stretches(1) + m_stretches(2)*m_stretches(2);
1893 T d2I_1 = 2*delta(a,b);
1894
1895 T d2I_2 = idelta(a,b)*4.0*m_stretches(a)*m_stretches(b) + delta(a,b)*2.0*(I_1 - m_stretches(a)*m_stretches(a));
1896
1897 T c2 = mu/(m_parvals.at(2)+1);
1898 T c1 = m_parvals.at(2)*c2;
1899
1900 return c1/2.0 * d2I_1 + c2/2.0 * d2I_2;
1901 }
1902
1903 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1904 template <enum Material _mat, bool _comp>
1905 typename std::enable_if<_comp && (_mat==Material::OG), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1906 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1907 {
1908 T tmp = 0.0;
1909 int n = (m_numPars-2)/2;
1910 T alpha_i, mu_i, Lambda;
1911 for (index_t k=0; k!=n; k++)
1912 {
1913 alpha_i = m_parvals.at(2*(k+1)+1);
1914 mu_i = m_parvals.at(2*(k+1));
1915 Lambda = math::pow(m_stretches(0),alpha_i) + math::pow(m_stretches(1),alpha_i) + math::pow(m_stretches(2),alpha_i);
1916 tmp += mu_i * math::pow(m_J_sq,-alpha_i/6.0) *
1917 ( - alpha_i/3. * ( math::pow(m_stretches(a),alpha_i-1.0) / m_stretches(b) + math::pow(m_stretches(b),alpha_i-1.0) / m_stretches(a)
1918 - 1./3. * 1. / (m_stretches(a)*m_stretches(b)) * Lambda )
1919 + delta(a,b) * ( (alpha_i - 1.) * math::pow(m_stretches(a),alpha_i-2.0) + Lambda / 3. * math::pow(m_stretches(a),-2.0) )
1920 );
1921 }
1922 tmp += _d2Psi_dab_vol(a,b);
1923 return tmp;
1924 }
1925
1926 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1927 template <enum Material _mat, bool _comp>
1928 typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1929 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1930 {
1931 T tmp = 0.0;
1932 int n = (m_numPars-2)/2;
1933 T alpha_i, mu_i;
1934 for (index_t k=0; k!=n; k++)
1935 {
1936 alpha_i = m_parvals.at(2*(k+1)+1);
1937 mu_i = m_parvals.at(2*(k+1));
1938 tmp += mu_i*math::pow(m_stretches(a),alpha_i-2)*(alpha_i-1)*delta(a,b);
1939 }
1940 return tmp;
1941 }
1942
1943 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1944 template <enum Material _mat, bool _comp>
1945 typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
_d2Psi_dab_impl(const index_t a,const index_t b) const1946 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
1947 {
1948 T mu = m_parvals.at(0) / (2 * (1 + m_parvals.at(1)));
1949 T d2I_1 = 2*delta(a,b);
1950 T lambda = m_parvals.at(0) * m_parvals.at(1) / ( (1. + m_parvals.at(1))*(1.-2.*m_parvals.at(1)));
1951
1952 return mu/2.0 * d2I_1 + mu * delta(a,b) / ( m_stretches(a) * m_stretches(b) ) + lambda / (2*m_stretches(a)*m_stretches(b)) * ( 2*m_J_sq - delta(a,b) * (m_J_sq - 1.0) );
1953 }
1954
1955 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1956 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi_dab_vol(const index_t a, const index_t b) const
1957 {
1958 GISMO_ENSURE(3 - 6 * m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1959 m_J_sq = math::pow(m_stretches(0)*m_stretches(1)*m_stretches(2),2.0);
1960 T beta = -2.0;
1961 T K = m_parvals.at(0) / ( 3 - 6 * m_parvals.at(1));
1962 return K / (beta*m_stretches(a)*m_stretches(b)) * ( beta*math::pow(m_J_sq,-beta/2.0) + delta(a,b) * (math::pow(m_J_sq,-beta/2.0) - 1.0) );
1963
1964 }
1965
1966 //--------------------------------------------------------------------------------------------------------------------------------------
1967
1968 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1969 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dJ_da(const index_t a) const
1970 {
1971 return 1.0/m_stretches(a);
1972 }
1973
1974 //--------------------------------------------------------------------------------------------------------------------------------------
1975
1976 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1977 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2J_dab(const index_t a, const index_t b) const
1978 {
1979 if (a==b)
1980 return 0.0;
1981 else
1982 return 1.0 / ( m_stretches(a) * m_stretches(b) );
1983 // return ( 1.0 - delta(a,b) ) / ( m_stretches(a) * m_stretches(b) );
1984 }
1985
1986 //--------------------------------------------------------------------------------------------------------------------------------------
1987
1988 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1989 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_p() const
1990 {
1991 return m_stretches(2) * _dPsi_da(2);
1992 }
1993
1994 //--------------------------------------------------------------------------------------------------------------------------------------
1995
1996 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1997 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dp_da(const index_t a) const
1998 {
1999 if (a==2)
2000 return m_stretches(2) * _d2Psi_dab(2,a) + _dPsi_da(2);
2001 else
2002 return m_stretches(2) * _d2Psi_dab(2,a);
2003
2004 // return m_stretches(2) * _d2Psi_dab(2,a) + delta(a,2) * _dPsi_da(2);
2005 }
2006
2007 //--------------------------------------------------------------------------------------------------------------------------------------
2008
2009 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2010 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sa(const index_t a) const
2011 {
2012 return _Sa_impl<comp>(a);
2013 }
2014
2015 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2016 template <bool _comp>
2017 typename std::enable_if<_comp, T>::type
_Sa_impl(const index_t a) const2018 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sa_impl(const index_t a) const
2019 {
2020 return 1.0/m_stretches(a) * _dPsi_da(a);
2021 }
2022
2023 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2024 template <bool _comp>
2025 typename std::enable_if<!_comp, T>::type
_Sa_impl(const index_t a) const2026 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Sa_impl(const index_t a) const
2027 {
2028 return 1.0/m_stretches(a) * (_dPsi_da(a) - _p() * _dJ_da(a) );
2029 }
2030
2031 //--------------------------------------------------------------------------------------------------------------------------------------
2032
2033 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2034 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dSa_db(const index_t a, const index_t b) const
2035 {
2036 return _dSa_db_impl<comp>(a,b);
2037 }
2038
2039 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2040 template <bool _comp>
2041 typename std::enable_if<_comp, T>::type
_dSa_db_impl(const index_t a,const index_t b) const2042 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dSa_db_impl(const index_t a, const index_t b) const
2043 {
2044 T tmp = 1.0/m_stretches(a) * _d2Psi_dab(a,b);
2045 if (a==b)
2046 tmp += - 1.0 / math::pow(m_stretches(a),2) * _dPsi_da(a);
2047 return tmp;
2048 }
2049
2050 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2051 template <bool _comp>
2052 typename std::enable_if<!_comp, T>::type
_dSa_db_impl(const index_t a,const index_t b) const2053 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dSa_db_impl(const index_t a, const index_t b) const
2054 {
2055 T tmp = 1.0/m_stretches(a) * ( _d2Psi_dab(a,b) - _dp_da(a)*_dJ_da(b) - _dp_da(b)*_dJ_da(a) - _p() * _d2J_dab(a,b) );
2056 if (a==b)
2057 tmp += - 1.0 / math::pow(m_stretches(a),2) * (_dPsi_da(a) - _p() * _dJ_da(a));
2058 return tmp;
2059 }
2060
2061 //--------------------------------------------------------------------------------------------------------------------------------------
2062
2063 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2064 T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cabcd(const index_t a, const index_t b, const index_t c, const index_t d) const
2065 {
2066 return _Cabcd_impl<comp>(a,b,c,d);
2067 }
2068
2069 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2070 template <bool _comp>
2071 typename std::enable_if<_comp, T>::type
_Cabcd_impl(const index_t a,const index_t b,const index_t c,const index_t d) const2072 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cabcd_impl(const index_t a, const index_t b, const index_t c, const index_t d) const
2073 {
2074 // Compute part with stress tensor involved.
2075 T frac = 0.0;
2076 T tmp = 0.0;
2077
2078 if (abs((m_stretches(a) - m_stretches(b)) / m_stretches(a)) < 1e-14)
2079 {
2080 // gsDebug<<"Stretches are equal; (abs((m_stretches(a) - m_stretches(b)) / m_stretches(a)) = "<<abs((m_stretches(a) - m_stretches(b)) / m_stretches(a))<<"\n";
2081 frac = 1.0 / (2.0 * m_stretches(a) ) * ( _dSa_db(b,b) - _dSa_db(a,b));
2082 }
2083 else
2084 frac = ( _Sa(b)-_Sa(a) ) / (math::pow(m_stretches(b),2) - math::pow(m_stretches(a),2));
2085
2086 GISMO_ENSURE( ( (a < 3) && (b < 3) && (c < 3) && (d < 3) ) , "Index out of range. a="<<a<<", b="<<b<<", c="<<c<<", d="<<d);
2087 if ( ( (a==b) && (c==d)) )
2088 tmp = 1/m_stretches(c) * _dSa_db(a,c);
2089 else if (( (a==d) && (b==c) && (a!=b) ) || ( ( (a==c) && (b==d) && (a!=b)) ))
2090 tmp = frac;
2091 // return 1/m_stretches(c) * _dSa_db(a,c) * delta(a,b) * delta(c,d) + frac * (delta(a,c)*delta(b,d) + delta(a,d)*delta(b,c)) * (1-delta(a,b));
2092
2093 return tmp;
2094 }
2095
2096 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2097 template <bool _comp>
2098 typename std::enable_if<!_comp, T>::type
_Cabcd_impl(const index_t a,const index_t b,const index_t c,const index_t d) const2099 gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_Cabcd_impl(const index_t a, const index_t b, const index_t c, const index_t d) const
2100 {
2101 // Compute part with stress tensor involved.
2102 T frac = 0.0;
2103 T tmp = 0.0;
2104
2105 if (abs((m_stretches(a) - m_stretches(b)) / m_stretches(a)) < 1e-14)
2106 {
2107 // gsDebug<<"Stretches are equal; (abs((m_stretches(a) - m_stretches(b)) / m_stretches(a)) = "<<abs((m_stretches(a) - m_stretches(b)) / m_stretches(a))<<"\n";
2108 frac = 1.0 / (2.0 * m_stretches(a) ) * ( _dSa_db(b,b) - _dSa_db(a,b));
2109 }
2110 else
2111 frac = ( _Sa(b)-_Sa(a) ) / (math::pow(m_stretches(b),2) - math::pow(m_stretches(a),2));
2112
2113 GISMO_ENSURE( ( (a < 2) && (b < 2) && (c < 2) && (d < 2) ) , "Index out of range. a="<<a<<", b="<<b<<", c="<<c<<", d="<<d);
2114 if ( ( (a==b) && (c==d)) )
2115 tmp = 1/m_stretches(c) * _dSa_db(a,c) + 1/(math::pow(m_stretches(a),2) * math::pow(m_stretches(c),2)) * ( math::pow(m_stretches(2),2) * _d2Psi_dab(2,2) + 2*_dPsi_da(2)*m_stretches(2) );
2116 else if (( (a==d) && (b==c) && (a!=b) ) || ( ( (a==c) && (b==d) && (a!=b)) ))
2117 tmp = frac;
2118 // return 1/m_stretches(c) * _dSa_db(a,c) * delta(a,b) * delta(c,d) + frac * (delta(a,c)*delta(b,d) + delta(a,d)*delta(b,c)) * (1-delta(a,b))
2119 // + delta(a,b)*delta(c,d)*1/(math::pow(m_stretches(a),2) * math::pow(m_stretches(c),2)) * ( math::pow(m_stretches(2),2) * _d2Psi_dab(2,2) + 2*_dPsi_da(2)*m_stretches(2) );
2120
2121 return tmp;
2122 }
2123
2124 //--------------------------------------------------------------------------------------------------------------------------------------
2125
2126 // template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2127 // T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_dPsi(const index_t a) const
2128 // {
2129 // T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
2130
2131 // GISMO_ENSURE( ( (a < 3) ) , "Index out of range. a="<<a);
2132 // GISMO_ENSURE(!comp,"Material model is not incompressible?");
2133
2134 // if (m_material==9)
2135 // {
2136 // return mu * m_stretches(a,0);
2137 // }
2138 // else if (m_material==3)
2139 // {
2140 // // return 2.0*m_parvals.at(0)*m_J0*m_G(i,j)*m_G(k,l) + m_G(i,k)*m_G(j,l) + m_G(i,l)*m_G(j,k);
2141 // }
2142 // else if (m_material==4)
2143 // {
2144 // gsMatrix<T> C(3,3);
2145 // C.setZero();
2146 // C.block(0,0,2,2) = m_Gcov_def.block(0,0,2,2);
2147 // C(2,2) = math::pow(m_J0,-2.0);
2148 // this->_computeStretch(C);
2149
2150 // return mu*m_stretches.at(a);
2151 // }
2152 // else
2153 // GISMO_ERROR("Material model not implemented.");
2154 // }
2155
2156 // template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
2157 // T gsMaterialMatrix<dim,T,matId,comp,mat,imp>::_d2Psi(const index_t a, const index_t b) const
2158 // {
2159 // T mu = m_parvals.at(0) / (2. * (1. + m_parvals.at(1)));
2160
2161 // GISMO_ENSURE( ( (a < 3) && (b < 3) ) , "Index out of range. a="<<a<<", b="<<b);
2162 // GISMO_ENSURE(!comp,"Material model is not incompressible?");
2163
2164 // if (m_material==9)
2165 // {
2166 // return ( a==b ? mu : 0.0 );
2167 // }
2168 // else if (m_material==3)
2169 // {
2170 // // return 2.0*m_parvals.at(0)*m_J0*m_G(i,j)*m_G(k,l) + m_G(i,k)*m_G(j,l) + m_G(i,l)*m_G(j,k);
2171 // }
2172 // else if (m_material==4)
2173 // {
2174 // if (a==b)
2175 // {
2176 // return mu;
2177 // }
2178 // else
2179 // return 0.0;
2180 // }
2181 // else
2182 // GISMO_ERROR("Material model not implemented.");
2183 // }
2184
2185 } // end namespace
2186