G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMaterialMatrixNonlinear.hpp
Go to the documentation of this file.
1
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#pragma once
23
24#include <gsKLShell/src/gsMaterialMatrixNonlinear.h>
25#include <gsKLShell/src/gsMaterialMatrixXml.hpp>
26#include <gsCore/gsFunction.h>
27
28using namespace gismo;
29
30template <short_t d, typename T>
31class Cfun : public gsFunction<T>
32{
33public:
34 Cfun(const gsMaterialMatrixBase<T> * materialMat, index_t patch, const gsVector<T> & u, const T & z)
35 :
36 m_materialMat(materialMat),
37 m_patchID(patch),
38 m_points(u),
39 m_z(z)
40 {
41
42 }
43
44 Cfun(const gsMaterialMatrixBase<T> * materialMat, index_t patch, const gsVector<T> & u)
45 :
46 m_materialMat(materialMat),
47 m_patchID(patch),
48 m_points(u)
49 {
50 m_z.resize(1,1);
51 m_z.setZero();
52 }
53
54 void eval_into(const gsMatrix<T>& C, gsMatrix<T>& result) const
55 {
56 result.resize(targetDim(),C.cols());
57 gsMatrix<T> Ctmp;
58 for (index_t k=0; k!=C.cols(); k++)
59 {
60 Ctmp = C.col(k);
61 Ctmp(2,0) *= 0.5;
62 result.col(k) = m_materialMat->eval3D_matrix_C(Ctmp,m_patchID,m_points,m_z,MaterialOutput::MatrixA);
63 }
64 }
65
66 short_t domainDim() const
67 {
68 return 3;
69 }
70
71 short_t targetDim() const
72 {
73 return 9;
74 }
75
76private:
77 const gsMaterialMatrixBase<T> * m_materialMat;
78 mutable index_t m_patchID;
79 mutable gsVector<T> m_points;
80 mutable T m_z;
81};
82
83namespace gismo
84{
85
86template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
87gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
88 const gsFunctionSet<T> & mp,
89 const gsFunctionSet<T> & thickness
90 )
91 :
92 Base(&mp,&thickness,nullptr)
93{
94 _initialize();
95}
96
97template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
98gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
99 const gsFunctionSet<T> & mp,
100 const gsFunctionSet<T> & thickness,
101 const std::vector<gsFunctionSet<T> *> &pars
102 )
103 :
104 gsMaterialMatrixNonlinear(&mp,&thickness,pars,nullptr)
105{}
106
107template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
108gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
109 const gsFunctionSet<T> & mp,
110 const gsFunctionSet<T> & thickness,
111 const std::vector<gsFunctionSet<T> *> &pars,
112 const gsFunctionSet<T> & Density
113 )
114 :
115 gsMaterialMatrixNonlinear(&mp,&thickness,pars,&Density)
116{}
117
118template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
119gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
120 const gsFunctionSet<T> & thickness,
121 const std::vector<gsFunctionSet<T> *> &pars
122 )
123 :
124 gsMaterialMatrixNonlinear(nullptr,&thickness,pars,nullptr)
125{}
126
127template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
128gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
129 const gsFunctionSet<T> & thickness,
130 const std::vector<gsFunctionSet<T> *> &pars,
131 const gsFunctionSet<T> & Density
132 )
133 :
134 gsMaterialMatrixNonlinear(nullptr,&thickness,pars,&Density)
135{}
136
137template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
138gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::gsMaterialMatrixNonlinear(
139 const gsFunctionSet<T> * mp,
140 const gsFunctionSet<T> * thickness,
141 const std::vector<gsFunctionSet<T> *> &pars,
142 const gsFunctionSet<T> * Density
143 )
144 :
145 Base(mp,thickness,Density)
146{
147 GISMO_ASSERT(pars.size()>=2,"Two or more material parameters should be assigned!");
148 this->setParameters(pars);
149 _initialize();
150}
151
152template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
153std::ostream & gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::print(std::ostream &os) const
154{
155 os <<"---------------------------------------------------------------------\n"
156 <<"---------------------Hyperelastic Material Info----------------------\n"
157 <<"---------------------------------------------------------------------\n\n";
158
159 os <<"Material model: \t";
160 if (comp)
161 os<<"Compressible ";
162 else
163 os<<"Incompressible ";
164
165 if (mat==Material::SvK)
166 os<<"Saint-Venant Kirchhoff";
167 else if (mat==Material::NH)
168 os<<"Neo-Hookean";
169 else if (mat==Material::MR)
170 os<<"Mooney-Rivlin";
171 else if (mat==Material::OG)
172 os<<"Ogden";
173 else if (mat==Material::NH_ext)
174 os<<"Neo-Hookean Extended";
175 else
176 gsWarn<<"Not specified";
177 os<<"\n";
178
179 os <<"Implementation: \t";
180 if (imp==Implementation::Analytical)
181 os<<"Analytical";
182 else if (imp==Implementation::Generalized)
183 os<<"Generalized";
184 else if (imp==Implementation::Spectral)
185 os<<"Spectral";
186 else
187 gsWarn<<"Not specified";
188 os<<" implementation\n";
189
190 os <<"---------------------------------------------------------------------\n\n";
191 return os;
192}
193
194template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
195void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::defaultOptions()
196{
197 Base::defaultOptions();
198 m_options.addInt("NumGauss","Number of Gaussian points through thickness",4);
199}
200
201template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
202void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_initialize()
203{
204 // Set default options
205 this->defaultOptions();
206}
207
208template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
209void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::pstretch_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
210{
211 this->_computePoints(patch,u);
212 _pstretch_into_impl<comp>(patch,u,result);
213}
214
215template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
216template <bool _comp>
217typename std::enable_if<!_comp, void>::type
218gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_pstretch_into_impl(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
219{
220 Base::pstretch_into(patch,u,result);
221}
222
223template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
224template <bool _comp>
225typename std::enable_if<_comp, void>::type
226gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_pstretch_into_impl(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
227{
228 result.resize(3, u.cols());
229 std::pair<gsVector<T>,gsMatrix<T>> res;
230
231 gsMatrix<T> zmat(1,1);
232 zmat<<0.0;
233 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,zmat);
234 T C33;
235 gsMatrix<T,3,3> c;
236 for (index_t i=0; i!= u.cols(); i++)
237 {
238 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
239 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,i);
240
241 this->_getMetric(i,0.0); // on point i, with height 0.0
242
243 C33 = C33s(0,i);
244
245 // Compute c
246 c.setZero();
247 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
248 c(2,2) = C33; // c33
249
250 res = this->_evalStretch(c,m_data.mine().m_gcon_ori);
251 result.col(i) = res.first;
252 }
253}
254
255
256template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
257void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::pstretchDir_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
258{
259 this->_computePoints(patch,u);
260 _pstretchDir_into_impl<comp>(patch,u,result);
261}
262
263template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
264template <bool _comp>
265typename std::enable_if<!_comp, void>::type
266gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_pstretchDir_into_impl(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
267{
268 Base::pstretchDir_into(patch,u,result);
269}
270
271template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
272template <bool _comp>
273typename std::enable_if<_comp, void>::type
274gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_pstretchDir_into_impl(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
275{
276 result.resize(9, u.cols());
277 std::pair<gsVector<T>,gsMatrix<T>> res;
278
279 gsMatrix<T> zmat(1,1);
280 zmat<<0.0;
281 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,zmat);
282 T C33;
283 gsMatrix<T,3,3> c;
284 for (index_t i=0; i!= u.cols(); i++)
285 {
286 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
287 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,i);
288
289 this->_getMetric(i,0.0); // on point i, with height 0.0
290
291 C33 = C33s(0,i);
292
293 // Compute c
294 c.setZero();
295 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
296 c(2,2) = C33; // c33
297
298 res = this->_evalStretch(c,m_data.mine().m_gcon_ori);
299 result.col(i) = res.second.reshape(9,1);
300 }
301}
302
303template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
304gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_matrix_C(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z, enum MaterialOutput /*out*/) const
305{
306 // Input: j index in-plane point
307 // z out-of-plane coordinate (through thickness) in R1 (z)
308 // Output: (n=u.cols(), m=z.rows())
309 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
310 return _eval3D_matrix_C_impl<mat,comp>(Cmat,patch,u,z);
311}
312
313template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
314template <enum Material _mat, bool _comp>
315typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
316gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
317{
318 this->_computePoints(patch,u);
319 gsMatrix<T> result = _eval3D_Incompressible_matrix_C(Cmat,patch, u, z);
320 return result;
321}
322
323template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
324template <enum Material _mat, bool _comp>
325typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
326gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
327{
328 this->_computePoints(patch,u);
329 gsMatrix<T> result = _eval3D_Incompressible_matrix_C(Cmat,patch, u, z);
330 return result;
331}
332
333template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
334template <enum Material _mat, bool _comp>
335typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
336gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
337{
338 this->_computePoints(patch,u);
339 gsMatrix<T> result = _eval3D_Compressible_matrix_C(Cmat,patch, u, z);
340 return result;
341}
342
343template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
344gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_dmatrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput /*out*/) const
345{
346 // return gsMatrix<T>::Zero(27,u.cols()*z.rows());
347
348 this->_computePoints(patch,u);
349
350 gsMatrix<T> result(27,u.cols()*z.rows());
351 result.setZero();
352 index_t colIdx;
353 for (index_t k=0; k!=u.cols(); k++)
354 {
355 // Evaluate material properties on the quadrature point
356 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
357 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
358
359 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
360 {
361 colIdx = j*u.cols()+k;
362 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
363 result.col(colIdx) = dCijkl(patch,u.col(k),z(j,k));
364 // gsAsMatrix<T,Dynamic,Dynamic> dCdC = result.reshapeCol(k,3,9);
365 // dCdC.row(2) *= 2;
366 }
367 }
368 return result;
369}
370
371template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
372inline gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl(const index_t patch, const gsVector<T> & u, const T z) const
373{
374 return dCijkl_impl<mat,comp>(patch,u,z);
375}
376
377template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
378template <enum Material _mat, bool _comp>
379inline typename std::enable_if<!(!_comp && _mat==Material::NH), gsMatrix<T>>::type
380gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl_impl(const index_t patch, const gsVector<T> & u, const T z) const
381{
382 gsMatrix<T> result;
383 Cfun<dim,T> Cfunc(this,patch,u,z);
384 gsMatrix<T> Cvoight(3,1);
385
386 Cvoight(0,0) = m_data.mine().m_Gcov_def(0,0);
387 Cvoight(1,0) = m_data.mine().m_Gcov_def(1,1);
388 Cvoight(2,0) = m_data.mine().m_Gcov_def(0,1);
389
390 Cfunc.deriv_into(Cvoight,result);
391 return result;
392}
393
394template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
395template <enum Material _mat, bool _comp>
396inline typename std::enable_if<(!_comp && _mat==Material::NH), gsMatrix<T>>::type
397gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl_impl(const index_t /*patch*/, const gsVector<T> & /*u*/, const T /*z*/) const
398{
399 gsMatrix<T> dCdC(3,9);
400 dCdC(0,0) = dCijkl_dCmn(0,0,0,0, 0,0); // C1111dC11
401 dCdC(1,0) = dCijkl_dCmn(0,0,0,0, 1,1); // C1111dC22
402 dCdC(2,0) = dCijkl_dCmn(0,0,0,0, 0,1); // C1111dC12
403
404 dCdC(0,1) = dCijkl_dCmn(0,0,1,1, 0,0); // C2211dC11
405 dCdC(1,1) = dCijkl_dCmn(0,0,1,1, 1,1); // C2211dC22
406 dCdC(2,1) = dCijkl_dCmn(0,0,1,1, 0,1); // C2211dC12
407
408 dCdC(0,2) = dCijkl_dCmn(0,0,0,1, 0,0); // C1211dC11
409 dCdC(1,2) = dCijkl_dCmn(0,0,0,1, 1,1); // C1211dC22
410 dCdC(2,2) = dCijkl_dCmn(0,0,0,1, 0,1); // C1211dC12
411
412 dCdC(0,3) = dCijkl_dCmn(1,1,0,0, 0,0); // C1122dC11
413 dCdC(1,3) = dCijkl_dCmn(1,1,0,0, 1,1); // C1122dC22
414 dCdC(2,3) = dCijkl_dCmn(1,1,0,0, 0,1); // C1122dC12
415
416 dCdC(0,4) = dCijkl_dCmn(1,1,1,1, 0,0); // C2222dC11
417 dCdC(1,4) = dCijkl_dCmn(1,1,1,1, 1,1); // C2222dC22
418 dCdC(2,4) = dCijkl_dCmn(1,1,1,1, 0,1); // C2222dC12
419
420 dCdC(0,5) = dCijkl_dCmn(1,1,0,1, 0,0); // C1222dC11
421 dCdC(1,5) = dCijkl_dCmn(1,1,0,1, 1,1); // C1222dC22
422 dCdC(2,5) = dCijkl_dCmn(1,1,0,1, 0,1); // C1222dC12
423
424 dCdC(0,6) = dCijkl_dCmn(0,1,0,0, 0,0); // C1112dC11
425 dCdC(1,6) = dCijkl_dCmn(0,1,0,0, 1,1); // C1112dC22
426 dCdC(2,6) = dCijkl_dCmn(0,1,0,0, 0,1); // C1112dC12
427
428 dCdC(0,7) = dCijkl_dCmn(0,1,1,1, 0,0); // C2212dC11
429 dCdC(1,7) = dCijkl_dCmn(0,1,1,1, 1,1); // C2212dC22
430 dCdC(2,7) = dCijkl_dCmn(0,1,1,1, 0,1); // C2212dC12
431
432 dCdC(0,8) = dCijkl_dCmn(0,1,0,1, 0,0); // C1212dC11
433 dCdC(1,8) = dCijkl_dCmn(0,1,0,1, 1,1); // C1212dC22
434 dCdC(2,8) = dCijkl_dCmn(0,1,0,1, 0,1); // C1212dC12
435 return dCdC.reshape(27,1);
436}
437
438template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
439inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl_dCmn(const index_t i, const index_t j, const index_t k, const index_t l, const index_t m, const index_t n) const
440{
441 return dCijkl_dCmn_impl<mat,comp>(i,j,k,l,m,n);
442}
443
444
445// template <class T>
446// T _dCinv_ij_dCmn(const gsMatrix<T> & Cinv, const index_t i, const index_t j, const index_t m, const index_t n)
447// {
448// return - 1/2 * ( Cinv(i,m) * Cinv(j,l) + Cinv(i,n) * Cinv(j,n) );
449// }
450
451template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
452inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dgconij_dCkl(const gsMatrix<T> & gcon, const index_t i, const index_t j, const index_t k, const index_t l) const
453{
454 return -1./2. * ( gcon(i,k) * gcon(j,l) + gcon(i,l) * gcon(j,k) );
455}
456
457template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
458template <enum Material _mat, bool _comp>
459inline typename std::enable_if<(!_comp && _mat==Material::NH), T>::type
460gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl_dCmn_impl(const index_t i, const index_t j, const index_t k, const index_t l, const index_t m, const index_t n) const
461{
462 // --------------------------
463 // Neo-Hookean
464 // --------------------------
465 const gsMatrix<T> Cinv = m_data.mine().m_Gcon_def;
466 const T J0_sq = m_data.mine().m_J0_sq;
467
468 const T X = (2.*Cinv(i,j)*Cinv(k,l) + Cinv(i,k)*Cinv(j,l) + Cinv(i,l)*Cinv(j,k));
469 const T dX= (2.*(_dgconij_dCkl(Cinv,i,j,m,n)*Cinv(k,l) + Cinv(i,j)*_dgconij_dCkl(Cinv,k,l,m,n))
470 + _dgconij_dCkl(Cinv,i,k,m,n)*Cinv(j,l) + Cinv(i,k)*_dgconij_dCkl(Cinv,j,l,m,n)
471 + _dgconij_dCkl(Cinv,i,l,m,n)*Cinv(j,k) + Cinv(i,l)*_dgconij_dCkl(Cinv,j,k,m,n)
472 );
473
474 const T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
475 return mu * (- 1. / J0_sq * Cinv(m,n) * X + 1. / J0_sq *dX);
476
477 // return -mu*1./m_data.mine().m_J0_sq*Cinv(m,n)*(2.*Cinv(i,j)*Cinv(k,l) + Cinv(i,k)*Cinv(j,l) + Cinv(i,l)*Cinv(j,k))
478 // +mu*1./m_data.mine().m_J0_sq*
479 // return mu*1./m_data.mine().m_J0_sq*(
480
481
482 // 2.*m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l)
483 // + m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l)
484 // + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k));
485
486
487 // gsMatrix<T> Cinv = m_data.mine().m_Gcon_def;
488 // return
489 // -mu*1./m_data.mine().m_J0_sq*Cinv(m,n)*(2.*Cinv(i,j)*Cinv(k,l) + Cinv(i,k)*Cinv(j,l) + Cinv(i,l)*Cinv(j,k))
490 // -mu *m_data.mine().m_J0_sq*1./2. *(
491 // (Cinv(j,k)*Cinv(l,n) + Cinv(j,l)*Cinv(k,n) + 2*Cinv(j,n)*Cinv(k,l))*Cinv(i,m)
492 // +(Cinv(j,k)*Cinv(l,m) + Cinv(j,l)*Cinv(k,m) + 2*Cinv(j,m)*Cinv(k,l))*Cinv(i,n)
493 // +(Cinv(i,k)*Cinv(l,n) + Cinv(i,l)*Cinv(k,n))*Cinv(j,m)
494 // +(Cinv(i,k)*Cinv(l,m) + Cinv(i,l)*Cinv(k,m))*Cinv(j,n)
495 // +2*(Cinv(k,m)*Cinv(l,n) + Cinv(k,n)*Cinv(l,m))*Cinv(i,j)
496 // );
497
498}
499
500template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
501template <enum Material _mat, bool _comp>
502inline typename std::enable_if<!(!_comp && _mat==Material::NH), T>::type
503gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::dCijkl_dCmn_impl(const index_t /*i*/, const index_t /*j*/, const index_t /*k*/, const index_t /*l*/, const index_t /*m*/, const index_t /*n*/) const
504{
505 GISMO_NO_IMPLEMENTATION;
506}
507
508template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
509gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput /*out*/) const
510{
511 // Input: j index in-plane point
512 // z out-of-plane coordinate (through thickness) in R1 (z)
513 // Output: (n=u.cols(), m=z.rows())
514 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
515 return _eval3D_matrix_impl<mat,comp>(patch,u,z);
516}
517
518template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
519template <enum Material _mat, bool _comp>
520typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
521gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
522{
523 this->_computePoints(patch,u);
524 gsMatrix<T> result = _eval3D_Incompressible_matrix(patch, u, z);
525 return result;
526}
527
528template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
529template <enum Material _mat, bool _comp>
530typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
531gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
532{
533 this->_computePoints(patch,u);
534 gsMatrix<T> result = _eval3D_Incompressible_matrix(patch, u, z);
535 return result;
536}
537
538template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
539template <enum Material _mat, bool _comp>
540typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
541gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
542{
543 this->_computePoints(patch,u);
544 gsMatrix<T> result = _eval3D_Compressible_matrix(patch, u, z);
545 return result;
546}
547
548template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
549gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
550{
551 return this->eval3D_stress(patch,u,z,out);
552}
553
554template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
555gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_vector_C(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z, enum MaterialOutput out) const
556{
557 return this->eval3D_stress_C(Cmat,patch,u,z,out);
558}
559
560template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
561gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_CauchyVector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
562{
563 return this->eval3D_CauchyStress(patch,u,z,out);
564}
565
566template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
567gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
568{
569 gsMatrix<T> Smat = eval3D_stress(patch,u,z,out);
570 gsMatrix<T> result(2, u.cols() * z.rows());
571 result.setZero();
572 gsMatrix<T,3,3> S;
573 index_t colIdx;
574 std::pair<gsVector<T>,gsMatrix<T>> res;
575 for (index_t k=0; k!=u.cols(); k++)
576 {
577 // Evaluate material properties on the quadrature point
578 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
579 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
580
581 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
582 {
583 colIdx = j * u.cols() + k;
584 S.setZero();
585 S(0,0) = Smat(0,colIdx);
586 S(1,1) = Smat(1,colIdx);
587 S(0,1) = S(1,0) = Smat(2,colIdx);
588 S(2,2) = 0;
589 res = this->_evalPStress(S);
590 result.col(colIdx) = res.first;
591 }
592 }
593 return result;
594}
595
596template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
597gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_pstressDir(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
598{
599 gsMatrix<T> Smat = eval3D_stress(patch,u,z,out);
600 gsMatrix<T> result(2, u.cols() * z.rows());
601 result.setZero();
602 gsMatrix<T,3,3> S;
603 index_t colIdx;
604 std::pair<gsVector<T>,gsMatrix<T>> res;
605 for (index_t k=0; k!=u.cols(); k++)
606 {
607 // Evaluate material properties on the quadrature point
608 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
609 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
610
611 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
612 {
613 colIdx = j * u.cols() + k;
614 S.setZero();
615 S(0,0) = Smat(0,colIdx);
616 S(1,1) = Smat(1,colIdx);
617 S(0,1) = S(1,0) = Smat(2,colIdx);
618 S(2,2) = 0;
619 res = this->_evalPStress(S);
620 result.col(colIdx) = res.second.reshape(9,1);
621 }
622 }
623 return result;
624}
625
626template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
627gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_CauchyPStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
628{
629 gsMatrix<T> Smat = eval3D_CauchyStress(patch,u,z,out);
630 gsMatrix<T> result(2, u.cols() * z.rows());
631 result.setZero();
632 gsMatrix<T,3,3> S;
633 index_t colIdx;
634 std::pair<gsVector<T>,gsMatrix<T>> res;
635 for (index_t k=0; k!=u.cols(); k++)
636 {
637 // Evaluate material properties on the quadrature point
638 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
639 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
640
641 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
642 {
643 colIdx = j * u.cols() + k;
644 S.setZero();
645 S(0,0) = Smat(0,colIdx);
646 S(1,1) = Smat(1,colIdx);
647 S(0,1) = S(1,0) = Smat(2,colIdx);
648 S(2,2) = 0;
649 res = this->_evalPStress(S);
650 result.col(j * u.cols() + k) = res.first;
651 }
652 }
653
654 return result;
655
656}
657
658template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
659gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_stress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput ) const
660{
661 this->_computePoints(patch,u);
662 return this->_eval3D_stress_impl<mat,comp>(patch,u,z);
663}
664
665template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
666template <enum Material _mat, bool _comp>
667typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
668gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
669{
670 gsMatrix<T> result = _eval3D_Incompressible_stress(patch, u, z);
671 return result;
672}
673
674template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
675template <enum Material _mat, bool _comp>
676typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
677gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
678{
679 gsMatrix<T> result = _eval3D_Incompressible_stress(patch, u, z);
680 return result;
681}
682
683template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
684template <enum Material _mat, bool _comp>
685typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
686gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
687{
688 gsMatrix<T> result = _eval3D_Compressible_stress(patch, u, z);
689 return result;
690}
691
692template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
693gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_stress_C(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z, enum MaterialOutput /*out*/) const
694{
695 this->_computePoints(patch,u);
696 return this->_eval3D_stress_C_impl<mat,comp>(Cmat,patch,u,z);
697}
698
699template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
700template <enum Material _mat, bool _comp>
701typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
702gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
703{
704 gsMatrix<T> result = _eval3D_Incompressible_stress_C(Cmat,patch, u, z);
705 return result;
706}
707
708template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
709template <enum Material _mat, bool _comp>
710typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
711gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
712{
713 gsMatrix<T> result = _eval3D_Incompressible_stress_C(Cmat,patch, u, z);
714 return result;
715}
716
717template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
718template <enum Material _mat, bool _comp>
719typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
720gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_stress_C_impl(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
721{
722 gsMatrix<T> result = _eval3D_Compressible_stress_C(Cmat,patch, u, z);
723 return result;
724}
725
726template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
727gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_stress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
728{
729 // Input: j index in-plane point
730 // z out-of-plane coordinate (through thickness) in R1 (z)
731 // Output: (n=u.cols(), m=z.rows())
732 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
733 gsMatrix<T> result(3, u.cols() * z.rows());
734 result.setZero();
735 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,z);
736 T C33;
737 gsMatrix<T,3,3> c, cinv;
738 index_t colIdx;
739 for (index_t k=0; k!=u.cols(); k++)
740 {
741 // Evaluate material properties on the quadrature point
742 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
743 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
744
745 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
746 {
747 colIdx = j*u.cols()+k;
748 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
749
750 // Avoid almost-zero entries for undeformed configurations
751 if (gsAllCloseRelativeToMax(m_data.mine().m_Gcon_ori,m_data.mine().m_Gcon_def,1e-12))
752 result.col(colIdx).setZero();
753 else
754 {
755 C33 = C33s(0,colIdx);
756
757 // Compute c
758 c.setZero();
759 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
760 c(2,2) = C33; // c33
761 cinv.setZero();
762 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
763 cinv(2,2) = 1.0/C33;
764 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * C33;
765
766 result(0,colIdx) = _Sij(0,0,c,cinv); // S11
767 result(1,colIdx) = _Sij(1,1,c,cinv); // S22
768 result(2,colIdx) = _Sij(0,1,c,cinv); // S12
769 }
770 }
771 }
772 return result;
773}
774
775template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
776gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_stress_C(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
777{
778 // Input: j index in-plane point
779 // z out-of-plane coordinate (through thickness) in R1 (z)
780 // Output: (n=u.cols(), m=z.rows())
781 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
782 gsMatrix<T> result(3, u.cols());
783 result.setZero();
784
785 gsMatrix<T> zmat(1,1);
786 zmat<<z;
787
788 gsMatrix<T> C33s = _eval3D_Compressible_C33(Cmat,patch,u,zmat);
789 T C33;
790 gsMatrix<T,3,3> c, cinv;
791
792 for (index_t k=0; k!=u.cols(); k++)
793 {
794 // Evaluate material properties on the quadrature point
795 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
796 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
797
798 this->_getMetric(0, z * m_data.mine().m_Tmat(0, 0), Cmat); // on point i, on height z(0,j)
799
800 // Avoid almost-zero entries for undeformed configurations
801 if (gsAllCloseRelativeToMax(m_data.mine().m_Gcon_ori,m_data.mine().m_Gcon_def,1e-12))
802 result.col(k).setZero();
803 else
804 {
805 C33 = C33s(0,k);
806
807 // Compute c
808 c.setZero();
809 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
810 c(2,2) = C33; // c33
811 cinv.setZero();
812 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
813 cinv(2,2) = 1.0/C33;
814 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * C33;
815
816 result(0,k) = _Sij(0,0,c,cinv); // S11
817 result(1,k) = _Sij(1,1,c,cinv); // S22
818 result(2,k) = _Sij(0,1,c,cinv); // S12
819 }
820 }
821 return result;
822}
823
824template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
825gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_stress(const index_t /*patch*/, const gsMatrix<T> & u, const gsMatrix<T>& z) const
826{
827 // gsInfo<<"TO DO: evaluate moments using thickness";
828 // Input: u in-plane points
829 // z matrix with, per point, a column with z integration points
830 // Output: (n=u.cols(), m=z.rows())
831 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
832 gsMatrix<T> result(3, u.cols() * z.rows());
833 result.setZero();
834 index_t colIdx;
835 for (index_t k=0; k!=u.cols(); k++)
836 {
837 // Evaluate material properties on the quadrature point
838 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
839 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
840
841 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
842 {
843 colIdx = j * u.cols() + k;
844 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
845
846 // Avoid almost-zero entries for undeformed configurations
847 if (gsAllCloseRelativeToMax(m_data.mine().m_Gcon_ori,m_data.mine().m_Gcon_def,1e-12))
848 result.col(colIdx).setZero();
849 else
850 {
851 result(0, colIdx) = _Sij(0, 0);
852 result(1, colIdx) = _Sij(1, 1);
853 result(2, colIdx) = _Sij(0, 1);
854 }
855 }
856 }
857 return result;
858}
859
860template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
861gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_stress_C(const gsMatrix<T> & Cmat, const index_t /*patch*/, const gsVector<T> & u, const T z) const
862{
863 // gsInfo<<"TO DO: evaluate moments using thickness";
864 // Input: u in-plane points
865 // z matrix with, per point, a column with z integration points
866 // Output: (n=u.cols(), m=z.rows())
867 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
868 gsMatrix<T> result(3, u.cols());
869 result.setZero();
870 for (index_t k=0; k!=u.cols(); k++)
871 {
872 // Evaluate material properties on the quadrature point
873 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
874 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
875
876 this->_getMetric(0, z * m_data.mine().m_Tmat(0, 0), Cmat); // on point i, on height z(0,j)
877
878 // Avoid almost-zero entries for undeformed configurations
879 if (gsAllCloseRelativeToMax(m_data.mine().m_Gcon_ori,m_data.mine().m_Gcon_def,1e-12))
880 result.col(k).setZero();
881 else
882 {
883 result(0, k) = _Sij(0, 0);
884 result(1, k) = _Sij(1, 1);
885 result(2, k) = _Sij(0, 1);
886 }
887 }
888 return result;
889}
890
891template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
892gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_detF(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput ) const
893{
894 this->_computePoints(patch,u);
895 return this->_eval3D_detF_impl<mat,comp>(patch,u,z);
896}
897
898template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
899template <enum Material _mat, bool _comp>
900typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
901gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_detF_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
902{
903 gsMatrix<T> result = _eval3D_Incompressible_detF(patch, u, z);
904 return result;
905}
906
907template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
908template <enum Material _mat, bool _comp>
909typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
910gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_detF_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
911{
912 gsMatrix<T> result = _eval3D_Incompressible_detF(patch, u, z);
913 return result;
914}
915
916template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
917template <enum Material _mat, bool _comp>
918typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
919gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_detF_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
920{
921 gsMatrix<T> result = _eval3D_Compressible_detF(patch, u, z);
922 return result;
923}
924
925template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
926gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_detF(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
927{
928 // Input: j index in-plane point
929 // z out-of-plane coordinate (through thickness) in R1 (z)
930 // Output: (n=u.cols(), m=z.rows())
931 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
932 gsMatrix<T> result(1, u.cols() * z.rows());
933 result.setZero();
934 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,z);
935 index_t colIdx;
936 for (index_t k=0; k!=u.cols(); k++)
937 {
938 // Evaluate material properties on the quadrature point
939 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
940 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
941
942 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
943 {
944 colIdx = j*u.cols()+k;
945 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
946 result(0,colIdx) = math::sqrt(m_data.mine().m_J0_sq*C33s(0,colIdx));
947 }
948 }
949 return result;
950}
951
952template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
953gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_detF(const index_t /*patch*/, const gsMatrix<T> & u, const gsMatrix<T>& z) const
954{
955 gsMatrix<T> result(1, u.cols() * z.rows());
956 result.setOnes();
957 return result;
958}
959
960template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
961gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::eval3D_CauchyStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput ) const
962{
963 this->_computePoints(patch,u);
964 return this->_eval3D_CauchyStress_impl<mat,comp>(patch,u,z);
965}
966
967template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
968template <enum Material _mat, bool _comp>
969typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
970gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_CauchyStress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
971{
972 gsMatrix<T> result = _eval3D_Incompressible_CauchyStress(patch, u, z);
973 return result;
974}
975
976template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
977template <enum Material _mat, bool _comp>
978typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
979gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_CauchyStress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
980{
981 gsMatrix<T> result = _eval3D_Incompressible_CauchyStress(patch, u, z);
982 return result;
983}
984
985template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
986template <enum Material _mat, bool _comp>
987typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
988gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_CauchyStress_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
989{
990 gsMatrix<T> result = _eval3D_Compressible_CauchyStress(patch, u, z);
991 return result;
992}
993
994template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
995gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_CauchyStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
996{
997 // Input: j index in-plane point
998 // z out-of-plane coordinate (through thickness) in R1 (z)
999 // Output: (n=u.cols(), m=z.rows())
1000 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1001
1002 gsMatrix<T> Smat = _eval3D_Compressible_stress(patch,u,z);
1003
1004 gsMatrix<T> result(3, u.cols() * z.rows());
1005 result.setZero();
1006 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,z);
1007 T C33;
1008 index_t colIdx;
1009 T detF;
1010 for (index_t k=0; k!=u.cols(); k++)
1011 {
1012 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1013 {
1014 colIdx = j * u.cols() + k;
1015 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
1016 C33 = C33s(0,colIdx);
1017 detF = math::sqrt(m_data.mine().m_J0_sq*C33);
1018 Smat.col(colIdx) /= detF;
1019 }
1020 }
1021 return result;
1022}
1023
1024template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1025gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_CauchyStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1026{
1027 // Same as _eval_Incompressible_vector, since J=1
1028 return _eval3D_Incompressible_stress(patch,u,z);
1029}
1030
1031template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1032void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::setYoungsModulus(const gsFunctionSet<T> & YoungsModulus)
1033{
1034 Base::setParameter(0,YoungsModulus);
1035}
1036
1037template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1038const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::getYoungsModulus() const
1039{
1040 return Base::getParameter(0);
1041}
1042
1043template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1044void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::setPoissonsRatio(const gsFunctionSet<T> & PoissonsRatio)
1045{
1046 Base::setParameter(1,PoissonsRatio);
1047}
1048
1049template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1050const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::getPoissonsRatio() const
1051{
1052 return Base::getParameter(1);
1053}
1054
1055template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1056template <enum Material _mat>
1057typename std::enable_if<_mat==Material::MR, void>::type
1058gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setRatio_impl(const gsFunctionSet<T> & Ratio)
1059{
1060 Base::setParameter(2,Ratio);
1061}
1062
1063template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1064template <enum Material _mat>
1065typename std::enable_if<_mat!=Material::MR, void>::type
1066gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setRatio_impl(const gsFunctionSet<T> & /*Ratio*/)
1067{
1068 GISMO_NO_IMPLEMENTATION;
1069}
1070
1071template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1072template <enum Material _mat>
1073typename std::enable_if<_mat==Material::MR, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1074gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getRatio_impl() const
1075{
1076 return Base::getParameter(2);
1077}
1078
1079template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1080template <enum Material _mat>
1081typename std::enable_if<_mat!=Material::MR, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1082gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getRatio_impl() const
1083{
1084 GISMO_NO_IMPLEMENTATION;
1085}
1086
1087template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1088template <enum Material _mat>
1089typename std::enable_if<_mat==Material::OG, void>::type
1090gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setMu_impl(const index_t & i, const gsFunctionSet<T> & Mu_i)
1091{
1092 Base::setParameter(2+2*i,Mu_i);
1093}
1094
1095template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1096template <enum Material _mat>
1097typename std::enable_if<_mat!=Material::OG, void>::type
1098gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setMu_impl(const index_t & /*i*/, const gsFunctionSet<T> & /*Mu_i*/)
1099{
1100 GISMO_NO_IMPLEMENTATION;
1101}
1102
1103template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1104template <enum Material _mat>
1105typename std::enable_if<_mat==Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1106gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getMu_impl(const index_t & i) const
1107{
1108 return Base::getParameter(2+2*i);
1109}
1110
1111template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1112template <enum Material _mat>
1113typename std::enable_if<_mat!=Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1114gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getMu_impl(const index_t & /*i*/) const
1115{
1116 GISMO_NO_IMPLEMENTATION;
1117}
1118
1119template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1120template <enum Material _mat>
1121typename std::enable_if<_mat==Material::OG, void>::type
1122gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setAlpha_impl(const index_t & i, const gsFunctionSet<T> & Alpha_i)
1123{
1124 Base::setParameter(2+2*i+1,Alpha_i);
1125}
1126
1127template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1128template <enum Material _mat>
1129typename std::enable_if<_mat!=Material::OG, void>::type
1130gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setAlpha_impl(const index_t & /*i*/, const gsFunctionSet<T> & /*Alpha_i*/)
1131{
1132 GISMO_NO_IMPLEMENTATION;
1133}
1134
1135template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1136template <enum Material _mat>
1137typename std::enable_if<_mat==Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1138gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getAlpha_impl(const index_t & i) const
1139{
1140 return Base::getParameter(2+2*i+1);
1141}
1142
1143template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1144template <enum Material _mat>
1145typename std::enable_if<_mat!=Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1146gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getAlpha_impl(const index_t & /*i*/) const
1147{
1148 GISMO_NO_IMPLEMENTATION;
1149}
1150
1151
1152template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1153gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_matrix_C(const gsMatrix<T> & Cmat, const index_t /*patch*/, const gsVector<T> & u, const T z) const
1154{
1155 // gsInfo<<"TO DO: evaluate moments using thickness";
1156 // Input: u in-plane points
1157 // z matrix with, per point, a column with z integration points
1158 // Output: (n=u.cols(), m=z.rows())
1159 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1160
1161 gsMatrix<T> result(9, u.cols());
1162 result.setZero();
1163
1164 for (index_t k=0; k!=u.cols(); k++)
1165 {
1166 // Evaluate material properties on the quadrature point
1167 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
1168 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
1169
1170 this->_getMetric(0, z * m_data.mine().m_Tmat(0, 0), Cmat); // on point i, on height z(0,j)
1171
1172 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(k,3,3);
1173 /*
1174 C = C1111, C1122, C1112
1175 symm, C2222, C2212
1176 symm, symm, C1212
1177 */
1178 C(0,0) = _Cijkl(0,0,0,0); // C1111
1179 C(1,1) = _Cijkl(1,1,1,1); // C2222
1180 C(2,2) = _Cijkl(0,1,0,1); // C1212
1181 C(1,0) = C(0,1) = _Cijkl(0,0,1,1); // C1122
1182 C(2,0) = C(0,2) = _Cijkl(0,0,0,1); // C1112
1183 C(2,1) = C(1,2) = _Cijkl(1,1,0,1); // C2212
1184 }
1185
1186 return result;
1187}
1188
1189template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1190gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Incompressible_matrix(const index_t /*patch*/, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1191{
1192 // gsInfo<<"TO DO: evaluate moments using thickness";
1193 // Input: u in-plane points
1194 // z matrix with, per point, a column with z integration points
1195 // Output: (n=u.cols(), m=z.rows())
1196 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1197
1198 gsMatrix<T> result(9, u.cols() * z.rows());
1199 result.setZero();
1200 index_t colIdx;
1201 for (index_t k=0; k!=u.cols(); k++)
1202 {
1203 // Evaluate material properties on the quadrature point
1204 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
1205 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
1206
1207 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1208 {
1209 colIdx = j*u.cols()+k;
1210 this->_getMetric(k,z(j,k) * m_data.mine().m_Tmat(0,k) ); // on point i, on height z(0,j)
1211
1212 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(colIdx,3,3);
1213 /*
1214 C = C1111, C1122, C1112
1215 symm, C2222, C2212
1216 symm, symm, C1212
1217 */
1218 C(0,0) = _Cijkl(0,0,0,0); // C1111
1219 C(1,1) = _Cijkl(1,1,1,1); // C2222
1220 C(2,2) = _Cijkl(0,1,0,1); // C1212
1221 C(1,0) = C(0,1) = _Cijkl(0,0,1,1); // C1122
1222 C(2,0) = C(0,2) = _Cijkl(0,0,0,1); // C1112
1223 C(2,1) = C(1,2) = _Cijkl(1,1,0,1); // C2212
1224 }
1225 }
1226
1227 return result;
1228}
1229
1230template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1231inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Cijkl(const index_t i, const index_t j, const index_t k, const index_t l) const
1232{
1233 GISMO_ENSURE( ( (i < 2) && (j < 2) && (k < 2) && (l < 2) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
1234 GISMO_ENSURE(!comp,"Material model is not incompressible?");
1235
1236 return _Cijkl_impl<mat,imp>(i,j,k,l);
1237}
1238
1239template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1240template <enum Material _mat, enum Implementation _imp>
1241inline typename std::enable_if<_mat==Material::SvK && _imp==Implementation::Analytical, T>::type
1242gsMaterialMatrixNonlinear<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
1243{
1244 // --------------------------
1245 // Saint Venant Kirchhoff
1246 // --------------------------
1247 T lambda, mu, Cconstant;
1248
1249 mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1250 GISMO_ENSURE((1.-2.*m_data.mine().m_parvals.at(1)) != 0, "Division by zero in construction of SvK material parameters! (1.-2.*m_data.mine().m_parvals.at(1)) = "<<(1.-2.*m_data.mine().m_parvals.at(1))<<"; m_data.mine().m_parvals.at(1) = "<<m_data.mine().m_parvals.at(1));
1251 lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1))) ;
1252 Cconstant = 2*lambda*mu/(lambda+2*mu);
1253
1254 return Cconstant*m_data.mine().m_Acon_ori(i,j)*m_data.mine().m_Acon_ori(k,l) + mu*(m_data.mine().m_Acon_ori(i,k)*m_data.mine().m_Acon_ori(j,l) + m_data.mine().m_Acon_ori(i,l)*m_data.mine().m_Acon_ori(j,k));
1255}
1256
1257template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1258template <enum Material _mat, enum Implementation _imp>
1259inline typename std::enable_if<_mat==Material::NH && _imp==Implementation::Analytical, T>::type
1260gsMaterialMatrixNonlinear<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
1261{
1262 // --------------------------
1263 // Neo-Hookean
1264 // --------------------------
1265 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1266 return mu*1./m_data.mine().m_J0_sq*(2.*m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l) + m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l) + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k));
1267}
1268
1269template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1270template <enum Material _mat, enum Implementation _imp>
1271inline typename std::enable_if<_mat==Material::MR && _imp==Implementation::Analytical, T>::type
1272gsMaterialMatrixNonlinear<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
1273{
1274 // --------------------------
1275 // Mooney-Rivlin
1276 // Parameter 3 is the ratio between c1 and c2.; c1 = m_data.mine().m_parvals.at(2)*c2
1277 // --------------------------
1278 GISMO_ENSURE(m_pars.size()==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1279 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1280 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1281 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1282 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1283
1284 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1285 T c2 = mu/(m_data.mine().m_parvals.at(2) + 1);
1286 T c1 = m_data.mine().m_parvals.at(2)*c2;
1287
1288 T Gabcd = - 1./2. * ( m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l) + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k) );
1289
1290 return (c1 + c2 * traceCt) *1./m_data.mine().m_J0_sq*(2.*m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l) + m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l) + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k))// correct
1291 - 2. * c2 / m_data.mine().m_J0_sq * ( m_data.mine().m_Gcon_ori(i,j) * m_data.mine().m_Gcon_def(k,l) + m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_ori(k,l)) // correct
1292 + 2. * c2 * m_data.mine().m_J0_sq * ( Gabcd + m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l) ); // Roohbakhshan
1293}
1294
1295template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1296template <enum Material _mat, enum Implementation _imp>
1297inline typename std::enable_if<_imp==Implementation::Spectral, T>::type
1298gsMaterialMatrixNonlinear<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
1299{
1300 // --------------------------
1301 // Stretch-based implementations
1302 // --------------------------
1303 T tmp = 0.0;
1304 gsMatrix<T> C(3,3);
1305 C.setZero();
1306 C.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1307 // C.block(0,0,2,2) = (m_gcov_def.transpose()*m_gcov_def).block(0,0,2,2);
1308 // gsDebugVar(m_gcov_def.transpose()*m_gcov_def);
1309 C(2,2) = 1./m_data.mine().m_J0_sq;
1310
1311 this->_computeStretch(C,m_data.mine().m_gcon_ori);
1312
1313 tmp = 0.0;
1314 for (index_t a = 0; a != 2; a++)
1315 {
1316 // C_iiii
1317 tmp += _Cabcd(a,a,a,a)*(
1318 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1319 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1320 );
1321
1322 for (index_t b = a+1; b != 2; b++)
1323 {
1324 // C_iijj = C_jjii
1325 tmp += _Cabcd(a,a,b,b)*(
1326 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1327 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1328 +
1329 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1330 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1331 );
1332
1333 // C_ijij = Cjiji = Cijji = Cjiij
1334 tmp += _Cabcd(a,b,a,b)*(
1335 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1336 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1337 +
1338 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1339 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1340 +
1341 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1342 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1343 +
1344 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1345 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1346 );
1347 }
1348 }
1349 return tmp;
1350}
1351
1352template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1353template <enum Material _mat, enum Implementation _imp>
1354inline typename std::enable_if<_imp==Implementation::Generalized, T>::type
1355gsMaterialMatrixNonlinear<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
1356{
1357 // --------------------------
1358 // General implementations
1359 // --------------------------
1360 return 4.0 * _d2Psi(i,j,k,l) + 4.0 * _d2Psi(2,2,2,2)*math::pow(m_data.mine().m_J0_sq,-2.0)*m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l)
1361 - 4.0/ m_data.mine().m_J0_sq * ( _d2Psi(2,2,i,j)*m_data.mine().m_Gcon_def(k,l) + _d2Psi(2,2,k,l)*m_data.mine().m_Gcon_def(i,j) )
1362 + 2.0 * _dPsi(2,2) / m_data.mine().m_J0_sq * (2.*m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l) + m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l) + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k));
1363}
1364
1365// template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1366// template <enum Material _mat>
1367// // OTHER CASES!
1368// typename std::enable_if<(_mat >= 30) && !(_mat==0) && !(_mat==2) && !(_mat==3), T>::type
1369// gsMaterialMatrixNonlinear<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
1370// {
1371// GISMO_ERROR("Material model unknown (model = "<<_mat<<"). Use gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::info() to see the options.");
1372// }
1373
1374 // else if (m_material==4)
1375 // GISMO_ERROR("Material model 4 is not invariant-based! Use 14 instead...");
1376 // else if (m_material==5)
1377 // GISMO_ERROR("Material model 5 is only for compressible materials...");
1378
1379
1380// Condensation of the 3D tensor for compressible materials
1381template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1382inline T gsMaterialMatrixNonlinear<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
1383{
1384 GISMO_ENSURE(c.cols()==c.rows(),"Matrix c must be square");
1385 GISMO_ENSURE(c.cols()==3,"Matrix c must be 3x3");
1386 GISMO_ENSURE(cinv.cols()==cinv.rows(),"Matrix cinv must be square");
1387 GISMO_ENSURE(cinv.cols()==3,"Matrix cinv must be 3x3");
1388 GISMO_ENSURE( ( (i <2) && (j <2) && (k <2) && (l <2) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
1389 GISMO_ENSURE(comp,"Material model is not compressible?");
1390
1391 return _Cijkl_impl<imp>(i,j,k,l,c,cinv);
1392}
1393
1394template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1395template <enum Implementation _imp>
1396inline typename std::enable_if< _imp==Implementation::Spectral, T>::type
1397gsMaterialMatrixNonlinear<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
1398{
1399 // static condensation is done before the projection
1400 this->_computeStretch(c,m_data.mine().m_gcon_ori);
1401 return _Cijkl3D(i,j,k,l,c,cinv);
1402}
1403
1404template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1405template <enum Implementation _imp>
1406inline typename std::enable_if<!(_imp==Implementation::Spectral), T>::type
1407gsMaterialMatrixNonlinear<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
1408{
1409 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);
1410}
1411
1412// 3D tensor for compressible materials
1413template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1414inline T gsMaterialMatrixNonlinear<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
1415{
1416 GISMO_ENSURE( ( (i <3) && (j <3) && (k <3) && (l <3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
1417 return _Cijkl3D_impl<mat,imp>(i,j,k,l,c,cinv);
1418}
1419
1420template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1421template <enum Material _mat, enum Implementation _imp>
1422inline typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
1423gsMaterialMatrixNonlinear<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
1424{
1425 GISMO_ERROR("Compressible material matrix requested, but not needed. How?");
1426}
1427
1428template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1429template <enum Material _mat, enum Implementation _imp>
1430inline typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1431gsMaterialMatrixNonlinear<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
1432{
1433 // --------------------------
1434 // Neo-Hookean
1435 // --------------------------
1436
1437 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
1438 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1439
1440 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
1441 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1442 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1443 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1444 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1445 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1446 T I_1 = traceCt + c(2,2);
1447 return 1.0 / 9.0 * mu * math::pow( m_data.mine().m_J_sq , -1.0/3.0 ) * ( 2.0 * I_1 * ( cinv(i,j)*cinv(k,l) - 3.0 * dCinv )
1448 - 6.0 *( m_data.mine().m_Gcon_ori(i,j)*cinv(k,l) + cinv(i,j)*m_data.mine().m_Gcon_ori(k,l) ) )
1449 + K * ( m_data.mine().m_J_sq*cinv(i,j)*cinv(k,l) + (m_data.mine().m_J_sq-1)*dCinv );
1450}
1451
1452template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1453template <enum Material _mat, enum Implementation _imp>
1454inline typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1455gsMaterialMatrixNonlinear<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
1456{
1457 // --------------------------
1458 // Mooney-Rivlin
1459 // --------------------------
1460 GISMO_ENSURE(m_pars.size()==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1461
1462 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
1463 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1464
1465 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
1466 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1467 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1468 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1469 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1470 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1471 T I_1 = traceCt + c(2,2);
1472 T I_2 = c(2,2) * traceCt + m_data.mine().m_J0_sq;
1473 T d2I_2 = idelta(i,2)*idelta(j,2)*idelta(k,2)*idelta(l,2)*( m_data.mine().m_J0_sq*( cinv(i,j)*cinv(k,l) + dCinv ) )
1474 + delta(i,2)*delta(j,2)*idelta(k,2)*idelta(l,2)*_dI_1(k,l)
1475 + idelta(i,2)*idelta(j,2)*delta(k,2)*delta(l,2)*_dI_1(i,j);
1476 T c2 = mu/(m_data.mine().m_parvals.at(2) + 1);
1477 T c1 = m_data.mine().m_parvals.at(2)*c2;
1478
1479 return 1.0/9.0 * c1 * math::pow(m_data.mine().m_J_sq, -1.0/3.0) * ( 2.0*I_1*cinv(i,j)*cinv(k,l) - 6.0*I_1*dCinv
1480 - 6.0*_dI_1(i,j)*cinv(k,l) - 6.0*cinv(i,j)*_dI_1(k,l) ) // + 9*d2I_1 = 0
1481 + 1.0/9.0 * c2 * math::pow(m_data.mine().m_J_sq, -2.0/3.0) * ( 8.0*I_2*cinv(i,j)*cinv(k,l) - 12.0*I_2*dCinv
1482 - 12.0*_dI_2(i,j,c,cinv)*cinv(k,l)- 12.0*cinv(i,j)*_dI_2(k,l,c,cinv)
1483 + 18.0*d2I_2 )
1484 + K * ( m_data.mine().m_J_sq*cinv(i,j)*cinv(k,l) + (m_data.mine().m_J_sq-1)*dCinv );
1485}
1486
1487template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1488template <enum Material _mat, enum Implementation _imp>
1489inline typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
1490gsMaterialMatrixNonlinear<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
1491{
1492 // --------------------------
1493 // Neo-Hookean 2
1494 // --------------------------
1495 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
1496 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
1497 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
1498
1499 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
1500 return - 2.0 * mu * dCinv + lambda * ( m_data.mine().m_J_sq*cinv(i,j)*cinv(k,l) + (m_data.mine().m_J_sq-1)*dCinv );
1501}
1502
1503template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1504template <enum Material _mat, enum Implementation _imp>
1505inline typename std::enable_if<_imp == Implementation::Spectral, T>::type
1506gsMaterialMatrixNonlinear<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
1507{
1508 // --------------------------
1509 // Stretch-based implementations
1510 // --------------------------
1511 T tmp;
1512 if ( (i==2) && (j==2) && (k==2) && (l==2) ) // if C3333 (no static condensation)
1513 tmp = _Cabcd(2,2,2,2);
1514 else
1515 {
1516 tmp = 0.0;
1517 T C = 0.0;
1518 T C2222 = _Cabcd(2,2,2,2);
1519 // T Cab22,C22ab;
1520 for (index_t a = 0; a != 2; a++)
1521 {
1522 // C_iiii
1523 // if (!((i==2 || j==2 || k==2 || l==2) && a!=2))
1524 // C = _Cabcd(a,a,a,a);
1525 C = _Cabcd(a,a,a,a) - math::pow(_Cabcd(2,2,a,a),2) / C2222;
1526 tmp += C*(
1527 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1528 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1529 );
1530
1531 for (index_t b = a+1; b != 2; b++)
1532 {
1533 // C_iijj
1534 // C = _Cabcd(a,a,b,b);
1535 C = _Cabcd(a,a,b,b) - _Cabcd(a,a,2,2) * _Cabcd(2,2,b,b) / C2222;
1536 tmp += C*(
1537 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1538 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1539 +
1540 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1541 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1542 );
1543
1544 // C_ijij = Cjiji = Cijji = Cjiij
1545 // C = _Cabcd(a,b,a,b);
1546 C = _Cabcd(a,b,a,b) - math::pow(_Cabcd(2,2,a,b),2) / C2222;
1547 tmp += C*(
1548 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1549 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1550 +
1551 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1552 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1553 +
1554 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(b)) )*
1555 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(a)) )
1556 +
1557 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(b)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )*
1558 ( m_data.mine().m_gcon_ori.col(k).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(l).dot(m_data.mine().m_stretchvec.col(b)) )
1559 );
1560 }
1561 }
1562 }
1563
1564 return tmp;
1565}
1566
1567template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1568template <enum Material _mat, enum Implementation _imp>
1569inline typename std::enable_if<_imp == Implementation::Generalized, T>::type
1570gsMaterialMatrixNonlinear<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
1571{
1572 // --------------------------
1573 // General implementations
1574 // --------------------------
1575 return 4.0 * _d2Psi(i,j,k,l,c,cinv);
1576}
1577
1578// Plan of approach:
1579// - Make function getMetricDeformed(k,z,basis,E), where k,z and basis are used to compute the undeformed metric
1580// - this new getMetricDeformed will compute the metric based on E, using m_Gcov_def etc =2E+m_Gcov_ori etc
1581// - Make an eval3D_matrix(u,z,...,E) which computes using E
1582// - Call that eval matrix inside the TFT given the u,z where the strain is also evaluated
1583
1584// template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
1585// gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_C(const gsMatrix<T> & Cmat) const
1586// {
1587// gsMatrix<T> result(9, u.cols() * z.rows());
1588// result.setZero();
1589
1590// for (index_t k=0; k!=u.cols(); k++)
1591// {
1592// // Evaluate material properties on the quadrature point
1593// for (index_t v=0; v!=m_parmat.rows(); v++)
1594// m_parvals.at(v) = m_parmat(v,k);
1595
1596// for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1597// {
1598// this->_getMetric(k, z(j,k) * m_data.mine().m_Tmat(0,k), Cmat);
1599
1600// }
1601// }
1602
1603// }
1604
1605/*
1606 Available class members:
1607 - m_data.mine().m_parvals.at(0)
1608 - m_data.mine().m_parvals.at(1)
1609 - m_metric
1610 - m_metric_def
1611 - m_data.mine().m_J0
1612 - m_data.mine().m_J
1613 - m_Cinv
1614*/
1615// template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1616// T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j) const { _Sij(i,j,NULL,NULL); }
1617
1618template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1619inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j) const
1620{
1621 return _Sij_impl<mat,imp>(i,j);
1622}
1623
1624template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1625template <enum Material _mat, enum Implementation _imp>
1626inline typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
1627gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1628{
1629 GISMO_NO_IMPLEMENTATION;
1630}
1631
1632template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1633template <enum Material _mat, enum Implementation _imp>
1634inline typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1635gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1636{
1637 // --------------------------
1638 // Neo-Hoookean
1639 // --------------------------
1640 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1641 return mu * (m_data.mine().m_Gcon_ori(i,j) - 1./m_data.mine().m_J0_sq * m_data.mine().m_Gcon_def(i,j) );
1642}
1643
1644template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1645template <enum Material _mat, enum Implementation _imp>
1646inline typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1647gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1648{
1649 // --------------------------
1650 // Mooney-Rivlin
1651 // Parameter 3 is the ratio between c1 and c2.
1652 // --------------------------
1653 GISMO_ENSURE(m_pars.size()==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1654 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1655 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1656 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1657 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1658 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1659
1660 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
1661 T c1 = m_data.mine().m_parvals.at(2)*c2;
1662
1663 return c1 * ( m_data.mine().m_Gcon_ori(i,j) - 1/m_data.mine().m_J0_sq * m_data.mine().m_Gcon_def(i,j) )
1664 + c2 / m_data.mine().m_J0_sq * (m_data.mine().m_Gcon_ori(i,j) - traceCt * m_data.mine().m_Gcon_def(i,j) ) + c2 * m_data.mine().m_J0_sq * m_data.mine().m_Gcon_def(i,j);
1665}
1666
1667template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1668template <enum Material _mat, enum Implementation _imp>
1669inline typename std::enable_if<_imp == Implementation::Spectral, T>::type
1670gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1671{
1672 T tmp = 0.0;
1673 gsMatrix<T> C(3,3);
1674 C.setZero();
1675 C.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1676 C(2,2) = 1./m_data.mine().m_J0_sq;
1677
1678 this->_computeStretch(C,m_data.mine().m_gcon_ori);
1679
1680 for (index_t a = 0; a != 2; a++)
1681 {
1682 tmp += _Sa(a)*(
1683 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )
1684 );
1685 }
1686 return tmp;
1687}
1688
1689template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1690template <enum Material _mat, enum Implementation _imp>
1691inline typename std::enable_if<_imp == Implementation::Generalized, T>::type
1692gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1693{
1694 // --------------------------
1695 // Generalized
1696 // --------------------------
1697 return 2.0 * _dPsi(i,j) - 2.0 * _dPsi(2,2) * math::pow(m_data.mine().m_J0_sq,-1.0)*m_data.mine().m_Gcon_def(i,j);
1698}
1699
1700//--------------------------------------------------------------------------------------------------------------------------------------
1701
1702template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1703inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
1704{
1705 return _Sij_impl<mat,imp>(i,j,c,cinv);
1706}
1707
1708template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1709template <enum Material _mat, enum Implementation _imp>
1710inline typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1711gsMaterialMatrixNonlinear<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
1712{
1713 // --------------------------
1714 // Neo-Hoookean
1715 // --------------------------
1716 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1717 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
1718 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1719 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1720 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1721 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1722 T I_1 = traceCt + c(2,2);
1723
1724 return mu * math::pow( m_data.mine().m_J_sq , -1.0/3.0 ) * ( m_data.mine().m_Gcon_ori(i,j) - 1.0/3.0 * I_1 * cinv(i,j) )
1725 + K * 0.5 * ( m_data.mine().m_J_sq - 1.0 ) * cinv(i,j);
1726}
1727
1728template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1729template <enum Material _mat, enum Implementation _imp>
1730inline typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1731gsMaterialMatrixNonlinear<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
1732{
1733 // --------------------------
1734 // Mooney-Rivlin
1735 // Parameter 3 is the ratio between c1 and c2.
1736 // --------------------------
1737 GISMO_ENSURE(m_pars.size()==3,"Mooney-Rivlin model needs to be a 3 parameter model");
1738 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1739 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
1740 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
1741 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
1742 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
1743 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
1744 T I_1 = traceCt + c(2,2);
1745 T I_2 = c(2,2) * traceCt + m_data.mine().m_J0_sq;
1746
1747 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
1748 T c1 = m_data.mine().m_parvals.at(2)*c2;
1749
1750 return c1 * math::pow( m_data.mine().m_J_sq , -1.0/3.0 ) * ( m_data.mine().m_Gcon_ori(i,j) - 1.0/3.0 * I_1 * cinv(i,j) )
1751 + c2 * math::pow( m_data.mine().m_J_sq , -2.0/3.0 ) * ( _dI_2(i,j,c,cinv)- 2.0/3.0 * I_2 * cinv(i,j) )
1752 + K * 0.5 * ( m_data.mine().m_J_sq - 1.0 ) * cinv(i,j);
1753}
1754
1755template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1756template <enum Material _mat, enum Implementation _imp>
1757inline typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
1758gsMaterialMatrixNonlinear<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
1759{
1760 // --------------------------
1761 // Neo-Hookean 2
1762 // --------------------------
1763 T mu = m_data.mine().m_parvals.at(0) / (2.*(1. + m_data.mine().m_parvals.at(1)));
1764 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
1765 return mu * m_data.mine().m_Gcon_ori(i,j) - mu * cinv(i,j) + lambda / 2.0 * ( m_data.mine().m_J_sq - 1 ) * cinv(i,j);
1766}
1767
1768template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1769template <enum Material _mat, enum Implementation _imp>
1770inline typename std::enable_if<_imp == Implementation::Spectral, T>::type
1771gsMaterialMatrixNonlinear<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
1772{
1773 T tmp = 0.0;
1774 this->_computeStretch(c,m_data.mine().m_gcon_ori);
1775 for (index_t a = 0; a != 3; a++)
1776 {
1777 tmp += _Sa(a)*(
1778 ( m_data.mine().m_gcon_ori.col(i).dot(m_data.mine().m_stretchvec.col(a)) )*( m_data.mine().m_gcon_ori.col(j).dot(m_data.mine().m_stretchvec.col(a)) )
1779 );
1780 }
1781 return tmp;
1782}
1783
1784template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1785template <enum Material _mat, enum Implementation _imp>
1786inline typename std::enable_if<_imp == Implementation::Generalized, T>::type
1787gsMaterialMatrixNonlinear<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
1788{
1789 // --------------------------
1790 // Generalized
1791 // --------------------------
1792 return 2.0 * _dPsi(i,j,c,cinv);
1793}
1794
1795//--------------------------------------------------------------------------------------------------------------------------------------
1796
1797template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1798inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sii(const index_t i) const // principle stresses
1799{
1800 gsMatrix<T> C(3,3);
1801 C.setZero();
1802 C.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1803 C(2,2) = 1./m_data.mine().m_J0_sq;
1804 return _Sii(i,C);
1805}
1806
1807template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1808inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sii(const index_t i, const gsMatrix<T> & c) const
1809{
1810 this->_computeStretch(c,m_data.mine().m_gcon_ori);
1811 return _Sa(i);
1812}
1813
1814template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1815gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_C33(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1816{
1817 // Input: j index in-plane point
1818 // z out-of-plane coordinate (through thickness) in R1 (z)
1819 // Output: (n=u.cols(), m=z.rows())
1820 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1821 this->_computePoints(patch,u);
1822 gsMatrix<T> result(1, u.cols() * z.rows());
1823 result.setZero();
1824 index_t colIdx;
1825 for (index_t k=0; k!=u.cols(); k++)
1826 {
1827 // Evaluate material properties on the quadrature point
1828 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
1829 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
1830
1831 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1832 {
1833 colIdx = j*u.cols()+k;
1834 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
1835
1836 // Define objects
1837 gsMatrix<T,3,3> c, cinv;
1838 T S33, C3333, dc33;
1839 // T S33_old;
1840 S33 = 0.0;
1841 dc33 = 0.0;
1842 C3333 = 1.0;
1843
1844 index_t itmax = 100;
1845 T tol = 1e-10;
1846
1847 // Initialize c
1848 c.setZero();
1849 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1850 c(2,2) = math::pow(m_data.mine().m_J0_sq,-1.0); // c33
1851 // c(2,2) = 1.0; // c33
1852 cinv.setZero();
1853 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
1854 cinv(2,2) = 1.0/c(2,2);
1855
1856 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2);
1857 S33 = _Sij(2,2,c,cinv);
1858 // S33_old = (S33 == 0.0) ? 1.0 : S33;
1859 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
1860
1861 dc33 = -2. * S33 / C3333;
1862 for (index_t it = 0; it < itmax; it++)
1863 {
1864 c(2,2) += dc33;
1865
1866 //GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
1867 cinv(2,2) = 1.0/c(2,2);
1868
1869 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2) ;
1870
1871 S33 = _Sij(2,2,c,cinv);
1872 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
1873
1874 dc33 = -2. * S33 / C3333;
1875 if (math::lessthan(math::abs(dc33),tol))
1876 {
1877 result(0,colIdx) = c(2,2);
1878 break;
1879 }
1880 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<", dc33 = "<<dc33<<" and tolerance = "<<tol<<"\n");
1881 }
1882 }
1883 }
1884 return result;
1885}
1886
1887
1888template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1889gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_C33(const gsMatrix<T> & Cmat, const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
1890{
1891 // Input: j index in-plane point
1892 // z out-of-plane coordinate (through thickness) in R1 (z)
1893 // Output: (n=u.cols(), m=z.rows())
1894 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1895 this->_computePoints(patch,u);
1896 gsMatrix<T> result(1, u.cols() * z.rows());
1897 result.setZero();
1898 index_t colIdx;
1899 for (index_t k=0; k!=u.cols(); k++)
1900 {
1901 // Evaluate material properties on the quadrature point
1902 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
1903 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
1904
1905 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
1906 {
1907 colIdx = j*u.cols()+k;
1908 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k),Cmat); // on point i, on height z(0,j)
1909
1910 // Define objects
1911 gsMatrix<T,3,3> c, cinv;
1912 T S33, C3333, dc33;
1913 // T S33_old;
1914 S33 = 0.0;
1915 dc33 = 0.0;
1916 C3333 = 1.0;
1917
1918 index_t itmax = 100;
1919 T tol = 1e-10;
1920
1921 // Initialize c
1922 c.setZero();
1923 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1924 c(2,2) = math::pow(m_data.mine().m_J0_sq,-1.0); // c33
1925 // c(2,2) = 1.0; // c33
1926 cinv.setZero();
1927 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
1928 cinv(2,2) = 1.0/c(2,2);
1929
1930 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2);
1931 S33 = _Sij(2,2,c,cinv);
1932 // S33_old = (S33 == 0.0) ? 1.0 : S33;
1933 C3333 = _Cijkl3D(2,2,2,2,c,cinv);
1934
1935 dc33 = -2. * S33 / C3333;
1936 for (index_t it = 0; it < itmax; it++)
1937 {
1938 c(2,2) += dc33;
1939
1940 //GISMO_ENSURE(c(2,2)>= 0,"ERROR in iteration "<<it<<"; c(2,2) = " << c(2,2) << " C3333=" << C3333 <<" S33=" << S33<<" dc33 = "<<dc33);
1941 cinv(2,2) = 1.0/c(2,2);
1942
1943 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * c(2,2) ;
1944
1945 S33 = _Sij(2,2,c,cinv);
1946 C3333 = _Cijkl3D(2,2,2,2,c,cinv); // or _Cijkl???
1947
1948 dc33 = -2. * S33 / C3333;
1949 if (math::lessthan(math::abs(dc33),tol))
1950 {
1951 result(0,colIdx) = c(2,2);
1952 break;
1953 }
1954 GISMO_ENSURE(it != itmax-1,"Error: Method did not converge, S33 = "<<S33<<", dc33 = "<<dc33<<" and tolerance = "<<tol<<"\n");
1955 }
1956 }
1957 }
1958 return result;
1959}
1960
1961template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1962gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_matrix_C(const gsMatrix<T> & Cmat, const index_t patch, const gsVector<T> & u, const T z) const
1963{
1964 // Input: j index in-plane point
1965 // z out-of-plane coordinate (through thickness) in R1 (z)
1966 // Output: (n=u.cols(), m=z.rows())
1967 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
1968 gsMatrix<T> result(9, u.cols());
1969 result.setZero();
1970 gsMatrix<T> zmat(1,1);
1971 zmat<<z;
1972 gsMatrix<T> C33s = _eval3D_Compressible_C33(Cmat,patch,u,zmat);
1973 T C33;
1974
1975 gsMatrix<T,3,3> c, cinv;
1976 for (index_t k=0; k!=u.cols(); k++)
1977 {
1978 // Evaluate material properties on the quadrature point
1979 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
1980 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
1981
1982 this->_getMetric(0, z * m_data.mine().m_Tmat(0, 0), Cmat); // on point i, on height z(0,j)
1983
1984 C33 = C33s(0,k);
1985
1986 // Compute c
1987 c.setZero();
1988 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
1989 c(2,2) = C33; // c33
1990 cinv.setZero();
1991 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
1992 cinv(2,2) = 1.0/C33;
1993 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * C33;
1994
1995 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(k,3,3);
1996 /*
1997 C = C1111, C1122, C1112
1998 symm, C2222, C2212
1999 symm, symm, C1212
2000 */
2001 C(0,0) = _Cijkl(0,0,0,0,c,cinv); // C1111
2002 C(1,1) = _Cijkl(1,1,1,1,c,cinv); // C2222
2003 C(2,2) = _Cijkl(0,1,0,1,c,cinv); // C1212
2004 C(1,0) = C(0,1) = _Cijkl(0,0,1,1,c,cinv); // C1122
2005 C(2,0) = C(0,2) = _Cijkl(0,0,0,1,c,cinv); // C1112
2006 C(2,1) = C(1,2) = _Cijkl(1,1,0,1,c,cinv); // C2212
2007 }
2008 return result;
2009}
2010
2011template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2012gsMatrix<T> gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_eval3D_Compressible_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
2013{
2014 // Input: j index in-plane point
2015 // z out-of-plane coordinate (through thickness) in R1 (z)
2016 // Output: (n=u.cols(), m=z.rows())
2017 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
2018 gsMatrix<T> result(9, u.cols() * z.rows());
2019 result.setZero();
2020 gsMatrix<T> C33s = _eval3D_Compressible_C33(patch,u,z);
2021 T C33;
2022 gsMatrix<T,3,3> c, cinv;
2023 index_t colIdx;
2024 for (index_t k=0; k!=u.cols(); k++)
2025 {
2026 // Evaluate material properties on the quadrature point
2027 for (index_t v=0; v!=m_data.mine().m_parmat.rows(); v++)
2028 m_data.mine().m_parvals.at(v) = m_data.mine().m_parmat(v,k);
2029
2030 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
2031 {
2032 colIdx = j*u.cols()+k;
2033 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
2034
2035 C33 = C33s(0,colIdx);
2036
2037 // Compute c
2038 c.setZero();
2039 c.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
2040 c(2,2) = C33; // c33
2041 cinv.setZero();
2042 cinv.block(0,0,2,2) = m_data.mine().m_Gcon_def.block(0,0,2,2);
2043 cinv(2,2) = 1.0/C33;
2044 m_data.mine().m_J_sq = m_data.mine().m_J0_sq * C33;
2045
2046 gsAsMatrix<T, Dynamic, Dynamic> C = result.reshapeCol(colIdx,3,3);
2047 /*
2048 C = C1111, C1122, C1112
2049 symm, C2222, C2212
2050 symm, symm, C1212
2051 */
2052 C(0,0) = _Cijkl(0,0,0,0,c,cinv); // C1111
2053 C(1,1) = _Cijkl(1,1,1,1,c,cinv); // C2222
2054 C(2,2) = _Cijkl(0,1,0,1,c,cinv); // C1212
2055 C(1,0) = C(0,1) = _Cijkl(0,0,1,1,c,cinv); // C1122
2056 C(2,0) = C(0,2) = _Cijkl(0,0,0,1,c,cinv); // C1112
2057 C(2,1) = C(1,2) = _Cijkl(1,1,0,1,c,cinv); // C2212
2058 }
2059 }
2060 return result;
2061}
2062
2063// ---------------------------------------------------------------------------------------------------------------------------------
2064// INCOMPRESSIBLE
2065// ---------------------------------------------------------------------------------------------------------------------------------
2066template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2067inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi(const index_t i, const index_t j) const
2068{
2069 GISMO_ENSURE( ( (i < 3) && (j < 3) ) , "Index out of range. i="<<i<<", j="<<j);
2070 GISMO_ENSURE(!comp,"Material model is not incompressible?");
2071 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
2072 return _dPsi_impl<mat>(i,j);
2073}
2074
2075template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2076template <enum Material _mat>
2077inline typename std::enable_if<_mat==Material::NH, T>::type
2078gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j) const
2079{
2080 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2081 return 0.5 * mu * m_data.mine().m_Gcon_ori(i,j);
2082}
2083
2084template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2085template <enum Material _mat>
2086inline typename std::enable_if<_mat==Material::MR, T>::type
2087gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_impl(const index_t i, const index_t j) const
2088{
2089 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2090 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2091 T c1 = m_data.mine().m_parvals.at(2)*c2;
2092 T tmp;
2093 if ((i==2) && (j==2))
2094 {
2095 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2096 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2097 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2098 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2099 tmp = c1/2.0 + c2 / 2.0 * traceCt;
2100 }
2101 else
2102 tmp = c1 / 2. * m_data.mine().m_Gcon_ori(i,j) + c2 / 2. * ( 1. / m_data.mine().m_J0_sq * m_data.mine().m_Gcon_ori(i,j) + m_data.mine().m_J0_sq * m_data.mine().m_Gcon_def(i,j) );
2103
2104 return tmp;
2105}
2106
2107
2108template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2109inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi(const index_t i, const index_t j, const index_t k, const index_t l) const
2110{
2111 GISMO_ENSURE( ( (i < 3) && (j < 3) && (k < 3) && (l < 3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
2112 GISMO_ENSURE(!comp,"Material model is not incompressible?");
2113 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
2114 return _d2Psi_impl<mat>(i,j,k,l);
2115}
2116
2117template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2118template <enum Material _mat>
2119inline typename std::enable_if<_mat==Material::NH, T>::type
2120gsMaterialMatrixNonlinear<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
2121{
2122 return 0.0;
2123}
2124
2125template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2126template <enum Material _mat>
2127inline typename std::enable_if<_mat==Material::MR, T>::type
2128gsMaterialMatrixNonlinear<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
2129{
2130 T tmp;
2131 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2132 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2133 // T c1 = m_data.mine().m_parvals.at(2)*c2;
2134 if ( ((i==2) && (j==2)) && !((k==2) || (l==2)) ) // _dPsi/d22dkl
2135 tmp = c2 / 2.0 * m_data.mine().m_Gcon_ori(k,l);
2136 else if ( !((i==2) && (j==2)) && ((k==2) || (l==2)) ) // _dPsi/dijd22
2137 tmp = c2 / 2.0 * m_data.mine().m_Gcon_ori(i,j);
2138 else if ( ((i==2) && (j==2)) && ((k==2) || (l==2)) ) // _dPsi/d22d22
2139 tmp = 0.0;
2140 else
2141 {
2142 T Gabcd = - 1./2. * ( m_data.mine().m_Gcon_def(i,k)*m_data.mine().m_Gcon_def(j,l) + m_data.mine().m_Gcon_def(i,l)*m_data.mine().m_Gcon_def(j,k) );
2143 tmp = c2 / 2.0 * m_data.mine().m_J0_sq * ( Gabcd + m_data.mine().m_Gcon_def(i,j)*m_data.mine().m_Gcon_def(k,l) );
2144 }
2145 return tmp;
2146}
2147
2148// ---------------------------------------------------------------------------------------------------------------------------------
2149// COMPRESSIBLE
2150// ---------------------------------------------------------------------------------------------------------------------------------
2151
2152//--------------------------------------------------------------------------------------------------------------------------------------
2153
2154template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2155inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dI_1(const index_t i, const index_t j) const
2156{
2157 return m_data.mine().m_Gcon_ori(i,j);
2158}
2159
2160//--------------------------------------------------------------------------------------------------------------------------------------
2161
2162template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2163inline T gsMaterialMatrixNonlinear<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
2164{
2165 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2166 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2167 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2168 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2169 return idelta(i,2)*idelta(j,2)*( c(2,2)*_dI_1(i,j) + m_data.mine().m_J0_sq*cinv(i,j) ) + delta(i,2)*delta(j,2)*traceCt;
2170}
2171
2172//--------------------------------------------------------------------------------------------------------------------------------------
2173
2174template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2175inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi(const index_t i, const index_t j, const gsMatrix<T> & c, const gsMatrix<T> & cinv) const
2176{
2177 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
2178 return _dPsi_impl<mat>(i,j,c,cinv);
2179}
2180
2181//--------------------------------------------------------------------------------------------------------------------------------------
2182
2183template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2184template <enum Material _mat>
2185inline typename std::enable_if<_mat==Material::NH, T>::type
2186gsMaterialMatrixNonlinear<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
2187{
2188 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2189 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2190 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2191 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2192 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2193 T I_1 = traceCt + c(2,2);
2194 return mu/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * ( - 1.0/3.0 * I_1 * cinv(i,j) + _dI_1(i,j) ) + _dPsi_vol(i,j,c,cinv);
2195}
2196
2197template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2198template <enum Material _mat>
2199inline typename std::enable_if<_mat==Material::MR, T>::type
2200gsMaterialMatrixNonlinear<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
2201{
2202 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2203 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2204 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2205 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2206 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2207 T I_1 = traceCt + c(2,2);
2208 T I_2 = c(2,2) * traceCt + m_data.mine().m_J0_sq;
2209 T c2= mu/(m_data.mine().m_parvals.at(2)+1);
2210 T c1= m_data.mine().m_parvals.at(2)*c2;
2211 return c1/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * ( - 1.0/3.0 * I_1 * cinv(i,j) + _dI_1(i,j) )
2212 + c2/2.0 * math::pow(m_data.mine().m_J_sq,-2./3.) * ( - 2.0/3.0 * I_2 * cinv(i,j) + _dI_2(i,j,c,cinv) )
2213 + _dPsi_vol(i,j,c,cinv);
2214}
2215
2216template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2217template <enum Material _mat>
2218inline typename std::enable_if<_mat==Material::NH_ext, T>::type
2219gsMaterialMatrixNonlinear<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
2220{
2221 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2222 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
2223 return mu / 2.0 * _dI_1(i,j) - mu / 2.0 * cinv(i,j) + lambda / 4.0 * ( m_data.mine().m_J_sq - 1 ) * cinv(i,j);
2224}
2225
2226// To do: add more models for volumetric part.
2227// Here, beta=2.
2228template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2229inline T gsMaterialMatrixNonlinear<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
2230{
2231 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
2232 return K * 0.25 * (m_data.mine().m_J_sq - 1.0) * cinv(i,j);
2233}
2234
2235//--------------------------------------------------------------------------------------------------------------------------------------
2236
2237template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2238inline T gsMaterialMatrixNonlinear<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
2239{
2240 GISMO_ENSURE( ( (i < 3) && (j < 3) && (k < 3) && (l < 3) ) , "Index out of range. i="<<i<<", j="<<j<<", k="<<k<<", l="<<l);
2241 GISMO_ENSURE(comp,"Material model is not compressible?");
2242 GISMO_ENSURE(imp==Implementation::Generalized,"Not generalized implementation");
2243 return _d2Psi_impl<mat>(i,j,k,l,c,cinv);
2244}
2245
2246template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2247template <enum Material _mat>
2248inline typename std::enable_if<_mat==Material::NH, T>::type
2249gsMaterialMatrixNonlinear<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
2250{
2251 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2252 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2253 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
2254 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2255 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2256 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2257 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2258 T I_1 = traceCt + c(2,2);
2259 return 1.0/9.0 * mu / 2.0 * math::pow(m_data.mine().m_J_sq, -1.0/3.0) * ( I_1*cinv(i,j)*cinv(k,l)
2260 - 3.0*_dI_1(i,j)*cinv(k,l) - 3.0*cinv(i,j)*_dI_1(k,l)
2261 - 3.0*I_1*dCinv )
2262 + _d2Psi_vol(i,j,k,l,c,cinv);
2263}
2264
2265template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2266template <enum Material _mat>
2267inline typename std::enable_if<_mat==Material::MR, T>::type
2268gsMaterialMatrixNonlinear<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
2269{
2270 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2271 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2272 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
2273 T traceCt = m_data.mine().m_Gcov_def(0,0)*m_data.mine().m_Gcon_ori(0,0) +
2274 m_data.mine().m_Gcov_def(0,1)*m_data.mine().m_Gcon_ori(0,1) +
2275 m_data.mine().m_Gcov_def(1,0)*m_data.mine().m_Gcon_ori(1,0) +
2276 m_data.mine().m_Gcov_def(1,1)*m_data.mine().m_Gcon_ori(1,1);
2277 T I_1 = traceCt + c(2,2);
2278 T I_2 = c(2,2) * traceCt + m_data.mine().m_J0_sq;
2279 T d2I_2 = idelta(i,2)*idelta(j,2)*idelta(k,2)*idelta(l,2)*( m_data.mine().m_J0_sq*( cinv(i,j)*cinv(k,l) + dCinv ) )
2280 + delta(i,2)*delta(j,2)*idelta(k,2)*idelta(l,2)*_dI_1(k,l)
2281 + idelta(i,2)*idelta(j,2)*delta(k,2)*delta(l,2)*_dI_1(i,j);
2282 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2283 T c1 = m_data.mine().m_parvals.at(2)*c2;
2284 // c1 = 0;
2285 return
2286 1.0/9.0 * c1 / 2.0 * math::pow(m_data.mine().m_J_sq, -1.0/3.0) * ( I_1*cinv(i,j)*cinv(k,l)
2287 - 3.0*_dI_1(i,j)*cinv(k,l) - 3.0*cinv(i,j)*_dI_1(k,l)
2288 - 3.0*I_1*dCinv ) // + 9*d2I_1 = 0
2289 + 1.0/9.0 * c2 / 2.0 * math::pow(m_data.mine().m_J_sq, -2.0/3.0) * ( 4.0*I_2*cinv(i,j)*cinv(k,l) - 6.0*I_2*dCinv
2290 - 6.0*_dI_2(i,j,c,cinv)*cinv(k,l)- 6.0*cinv(i,j)*_dI_2(k,l,c,cinv)
2291 + 9.0*d2I_2 )
2292 + _d2Psi_vol(i,j,k,l,c,cinv);
2293}
2294
2295template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2296template <enum Material _mat>
2297inline typename std::enable_if<_mat==Material::NH_ext, T>::type
2298gsMaterialMatrixNonlinear<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
2299{
2300 T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2301 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
2302 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
2303 return - mu / 2.0 * dCinv + lambda / 4.0 * ( m_data.mine().m_J_sq*cinv(i,j)*cinv(k,l) + (m_data.mine().m_J_sq-1.0)*dCinv );
2304}
2305
2306template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2307inline T gsMaterialMatrixNonlinear<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
2308{
2309 T dCinv = - 1./2.*( cinv(i,k)*cinv(j,l) + cinv(i,l)*cinv(j,k) );
2310 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
2311 return K * 0.25 * ( m_data.mine().m_J_sq*cinv(i,j)*cinv(k,l) + (m_data.mine().m_J_sq-1.0)*dCinv );
2312}
2313
2314// ---------------------------------------------------------------------------------------------------------------------------------
2315// STRETCHES
2316// ---------------------------------------------------------------------------------------------------------------------------------
2317
2318//--------------------------------------------------------------------------------------------------------------------------------------
2319
2320template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2321inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da(const index_t a) const
2322{
2323 GISMO_ENSURE( a < 3 , "Index out of range. a="<<a);
2324 return _dPsi_da_impl<mat,comp>(a);
2325}
2326
2327template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2328template <enum Material _mat, bool _comp>
2329inline typename std::enable_if<_comp && (_mat==Material::NH), T>::type
2330gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2331{
2332 GISMO_ENSURE( a < 3 , "Index out of range. a="<<a);
2333 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2334 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2335 T _dI_1a = 2*m_data.mine().m_stretches(a);
2336
2337 return mu/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * ( -2./3. * I_1 / m_data.mine().m_stretches(a) + _dI_1a )
2338 + _dPsi_da_vol(a);
2339}
2340
2341template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2342template <enum Material _mat, bool _comp>
2343inline typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
2344gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2345{
2346 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2347 T _dI_1a = 2*m_data.mine().m_stretches(a);
2348 return mu/2 * _dI_1a;
2349}
2350
2351template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2352template <enum Material _mat, bool _comp>
2353inline typename std::enable_if<_comp && (_mat==Material::MR), T>::type
2354gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2355{
2356 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2357 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2358 T _dI_1a = 2*m_data.mine().m_stretches(a);
2359 T I_2 = math::pow(m_data.mine().m_stretches(0),2.)*math::pow(m_data.mine().m_stretches(1),2.)
2360 + math::pow(m_data.mine().m_stretches(1),2.)*math::pow(m_data.mine().m_stretches(2),2.)
2361 + math::pow(m_data.mine().m_stretches(0),2.)*math::pow(m_data.mine().m_stretches(2),2.);
2362 T _dI_2a = 2*m_data.mine().m_stretches(a)*( I_1 - math::pow(m_data.mine().m_stretches(a),2.0) );
2363
2364 T c2= mu/(m_data.mine().m_parvals.at(2)+1);
2365 T c1= m_data.mine().m_parvals.at(2)*c2;
2366 return c1/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * ( -2./3. * I_1 / m_data.mine().m_stretches(a) + _dI_1a )
2367 + c2/2.0 * math::pow(m_data.mine().m_J_sq,-2./3.) * ( -4./3. * I_2 / m_data.mine().m_stretches(a) + _dI_2a )
2368 + _dPsi_da_vol(a);
2369}
2370
2371template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2372template <enum Material _mat, bool _comp>
2373inline typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
2374gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2375{
2376 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2377 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2378 T _dI_1a = 2*m_data.mine().m_stretches(a);
2379 T _dI_2a = 2*m_data.mine().m_stretches(a)*( I_1 - math::pow(m_data.mine().m_stretches(a),2.0) );
2380
2381 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2382 T c1 = m_data.mine().m_parvals.at(2)*c2;
2383 return c1/2.0*_dI_1a + c2/2.0*_dI_2a;
2384}
2385
2386template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2387template <enum Material _mat, bool _comp>
2388inline typename std::enable_if<_comp && (_mat==Material::OG), T>::type
2389gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2390{
2391 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2392 T tmp = 0.0;
2393 index_t n = (m_pars.size()-2)/2;
2394 T alpha_i, mu_i, Lambda;
2395 for (index_t k=0; k!=n; k++)
2396 {
2397 alpha_i = m_data.mine().m_parvals.at(2*(k+1)+1);
2398 mu_i = m_data.mine().m_parvals.at(2*(k+1));
2399 Lambda = math::pow(m_data.mine().m_stretches(0),alpha_i) + math::pow(m_data.mine().m_stretches(1),alpha_i) + math::pow(m_data.mine().m_stretches(2),alpha_i);
2400 tmp += mu_i * math::pow(m_data.mine().m_J_sq,-alpha_i/6.0) * ( math::pow(m_data.mine().m_stretches(a),alpha_i-1) - 1./3. * 1./m_data.mine().m_stretches(a) * Lambda );
2401 }
2402 return tmp + _dPsi_da_vol(a);
2403}
2404template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2405template <enum Material _mat, bool _comp>
2406inline typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
2407gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2408{
2409 T tmp = 0.0;
2410 index_t n = (m_pars.size()-2)/2;
2411 T alpha_i, mu_i;
2412 for (index_t k=0; k!=n; k++)
2413 {
2414 alpha_i = m_data.mine().m_parvals.at(2*(k+1)+1);
2415 mu_i = m_data.mine().m_parvals.at(2*(k+1));
2416 tmp += mu_i*math::pow(m_data.mine().m_stretches(a),alpha_i-1);
2417 }
2418 return tmp;
2419}
2420
2421template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2422template <enum Material _mat, bool _comp>
2423inline typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
2424gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_impl(const index_t a) const
2425{
2426 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2427 T _dI_1a = 2*m_data.mine().m_stretches(a);
2428 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2429 // choose compressibility function (and parameter)
2430 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
2431
2432 return mu/2.0 * _dI_1a - mu / m_data.mine().m_stretches(a) + lambda / (m_data.mine().m_stretches(a)*2) * (m_data.mine().m_J_sq-1.0);
2433}
2434
2435template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2436inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi_da_vol(const index_t a) const
2437{
2438 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2439 T beta = -2.0;
2440 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
2441 return K / (m_data.mine().m_stretches(a)*beta) * (1.0 - math::pow(m_data.mine().m_J_sq,-beta/2.0));
2442}
2443
2444//--------------------------------------------------------------------------------------------------------------------------------------
2445
2446template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2447inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab(const index_t a, const index_t b) const
2448{
2449 GISMO_ENSURE( ( (a < 3) && (b < 3) ) , "Index out of range. a="<<a<<", b="<<b);
2450 return _d2Psi_dab_impl<mat,comp>(a,b);
2451}
2452
2453template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2454template <enum Material _mat, bool _comp>
2455inline typename std::enable_if<_comp && (_mat==Material::NH), T>::type
2456gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2457{
2458 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2459
2460 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2461 T _dI_1a = 2*m_data.mine().m_stretches(a);
2462 T _dI_1b = 2*m_data.mine().m_stretches(b);
2463 T d2I_1 = 2*delta(a,b);
2464 return mu/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * (
2465 -2./3. * 1. / m_data.mine().m_stretches(b) * ( -2./3. * I_1 / m_data.mine().m_stretches(a) + _dI_1a )
2466 -2./3. * 1. / m_data.mine().m_stretches(a) * _dI_1b
2467 +d2I_1
2468 +2./3. * delta(a,b) * I_1 / (m_data.mine().m_stretches(a)*m_data.mine().m_stretches(a))
2469 )
2470 + _d2Psi_dab_vol(a,b);
2471}
2472
2473template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2474template <enum Material _mat, bool _comp>
2475inline typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
2476gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2477{
2478 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2479 T d2I_1 = 2*delta(a,b);
2480 return mu/2 * d2I_1;
2481}
2482
2483template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2484template <enum Material _mat, bool _comp>
2485inline typename std::enable_if<_comp && (_mat==Material::MR), T>::type
2486gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2487{
2488 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2489
2490 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2491 T _dI_1a = 2*m_data.mine().m_stretches(a);
2492 T _dI_1b = 2*m_data.mine().m_stretches(b);
2493 T d2I_1 = 2*delta(a,b);
2494
2495 T I_2 = math::pow(m_data.mine().m_stretches(0),2.)*math::pow(m_data.mine().m_stretches(1),2.)
2496 + math::pow(m_data.mine().m_stretches(1),2.)*math::pow(m_data.mine().m_stretches(2),2.)
2497 + math::pow(m_data.mine().m_stretches(0),2.)*math::pow(m_data.mine().m_stretches(2),2.);
2498 T _dI_2a = 2*m_data.mine().m_stretches(a)*( I_1 - math::pow(m_data.mine().m_stretches(a),2.0) );
2499 T _dI_2b = 2*m_data.mine().m_stretches(b)*( I_1 - math::pow(m_data.mine().m_stretches(b),2.0) );
2500 T d2I_2 = idelta(a,b)*4.0*m_data.mine().m_stretches(a)*m_data.mine().m_stretches(b) + delta(a,b)*2.0*(I_1 - m_data.mine().m_stretches(a)*m_data.mine().m_stretches(a));
2501
2502 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2503 T c1 = m_data.mine().m_parvals.at(2)*c2;
2504 return
2505 c1/2.0 * math::pow(m_data.mine().m_J_sq,-1./3.) * (
2506 -2./3. * 1. / m_data.mine().m_stretches(b) * ( -2./3. * I_1 / m_data.mine().m_stretches(a) + _dI_1a )
2507 -2./3. * 1. / m_data.mine().m_stretches(a) * _dI_1b
2508 +d2I_1
2509 +2./3. * delta(a,b) * I_1 / (m_data.mine().m_stretches(a)*m_data.mine().m_stretches(a))
2510 )
2511 + c2/2.0 * math::pow(m_data.mine().m_J_sq,-2./3.) * (
2512 -4./3. * 1. / m_data.mine().m_stretches(b) * ( -4./3. * I_2 / m_data.mine().m_stretches(a) + _dI_2a )
2513 -4./3. * 1. / m_data.mine().m_stretches(a) * _dI_2b
2514 +d2I_2
2515 +4./3. * delta(a,b) * I_2 / (m_data.mine().m_stretches(a)*m_data.mine().m_stretches(a))
2516 )
2517 + _d2Psi_dab_vol(a,b);
2518}
2519
2520template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2521template <enum Material _mat, bool _comp>
2522inline typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
2523gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2524{
2525 GISMO_ENSURE(m_pars.size()==3,"Mooney-Rivlin model needs to be a 3 parameter model");
2526 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2527
2528 T I_1 = m_data.mine().m_stretches(0)*m_data.mine().m_stretches(0) + m_data.mine().m_stretches(1)*m_data.mine().m_stretches(1) + m_data.mine().m_stretches(2)*m_data.mine().m_stretches(2);
2529 T d2I_1 = 2*delta(a,b);
2530
2531 T d2I_2 = idelta(a,b)*4.0*m_data.mine().m_stretches(a)*m_data.mine().m_stretches(b) + delta(a,b)*2.0*(I_1 - m_data.mine().m_stretches(a)*m_data.mine().m_stretches(a));
2532
2533 T c2 = mu/(m_data.mine().m_parvals.at(2)+1);
2534 T c1 = m_data.mine().m_parvals.at(2)*c2;
2535
2536 return c1/2.0 * d2I_1 + c2/2.0 * d2I_2;
2537}
2538
2539template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2540template <enum Material _mat, bool _comp>
2541inline typename std::enable_if<_comp && (_mat==Material::OG), T>::type
2542gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2543{
2544 T tmp = 0.0;
2545 index_t n = (m_pars.size()-2)/2;
2546 T alpha_i, mu_i, Lambda;
2547 for (index_t k=0; k!=n; k++)
2548 {
2549 alpha_i = m_data.mine().m_parvals.at(2*(k+1)+1);
2550 mu_i = m_data.mine().m_parvals.at(2*(k+1));
2551 Lambda = math::pow(m_data.mine().m_stretches(0),alpha_i) + math::pow(m_data.mine().m_stretches(1),alpha_i) + math::pow(m_data.mine().m_stretches(2),alpha_i);
2552 tmp += mu_i * math::pow(m_data.mine().m_J_sq,-alpha_i/6.0) *
2553 ( - alpha_i/3. * ( math::pow(m_data.mine().m_stretches(a),alpha_i-1.0) / m_data.mine().m_stretches(b) + math::pow(m_data.mine().m_stretches(b),alpha_i-1.0) / m_data.mine().m_stretches(a)
2554 - 1./3. * 1. / (m_data.mine().m_stretches(a)*m_data.mine().m_stretches(b)) * Lambda )
2555 + delta(a,b) * ( (alpha_i - 1.) * math::pow(m_data.mine().m_stretches(a),alpha_i-2.0) + Lambda / 3. * math::pow(m_data.mine().m_stretches(a),-2.0) )
2556 );
2557 }
2558 tmp += _d2Psi_dab_vol(a,b);
2559 return tmp;
2560}
2561
2562template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2563template <enum Material _mat, bool _comp>
2564inline typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
2565gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2566{
2567 T tmp = 0.0;
2568 index_t n = (m_pars.size()-2)/2;
2569 T alpha_i, mu_i;
2570 for (index_t k=0; k!=n; k++)
2571 {
2572 alpha_i = m_data.mine().m_parvals.at(2*(k+1)+1);
2573 mu_i = m_data.mine().m_parvals.at(2*(k+1));
2574 tmp += mu_i*math::pow(m_data.mine().m_stretches(a),alpha_i-2)*(alpha_i-1)*delta(a,b);
2575 }
2576 return tmp;
2577}
2578
2579template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2580template <enum Material _mat, bool _comp>
2581inline typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
2582gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_impl(const index_t a, const index_t b) const
2583{
2584 T mu = m_data.mine().m_parvals.at(0) / (2 * (1 + m_data.mine().m_parvals.at(1)));
2585 T d2I_1 = 2*delta(a,b);
2586 T lambda = m_data.mine().m_parvals.at(0) * m_data.mine().m_parvals.at(1) / ( (1. + m_data.mine().m_parvals.at(1))*(1.-2.*m_data.mine().m_parvals.at(1)));
2587
2588 return mu/2.0 * d2I_1 + mu * delta(a,b) / ( m_data.mine().m_stretches(a) * m_data.mine().m_stretches(b) ) + lambda / (2*m_data.mine().m_stretches(a)*m_data.mine().m_stretches(b)) * ( 2*m_data.mine().m_J_sq - delta(a,b) * (m_data.mine().m_J_sq - 1.0) );
2589}
2590
2591template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2592inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi_dab_vol(const index_t a, const index_t b) const
2593{
2594 GISMO_ENSURE(3 - 6 * m_data.mine().m_parvals.at(1) != 0, "Bulk modulus is infinity for compressible material model. Try to use incompressible models.");
2595 m_data.mine().m_J_sq = math::pow(m_data.mine().m_stretches(0)*m_data.mine().m_stretches(1)*m_data.mine().m_stretches(2),2.0);
2596 T beta = -2.0;
2597 T K = m_data.mine().m_parvals.at(0) / ( 3 - 6 * m_data.mine().m_parvals.at(1));
2598 return K / (beta*m_data.mine().m_stretches(a)*m_data.mine().m_stretches(b)) * ( beta*math::pow(m_data.mine().m_J_sq,-beta/2.0) + delta(a,b) * (math::pow(m_data.mine().m_J_sq,-beta/2.0) - 1.0) );
2599
2600}
2601
2602//--------------------------------------------------------------------------------------------------------------------------------------
2603
2604template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2605inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dJ_da(const index_t a) const
2606{
2607 return 1.0/m_data.mine().m_stretches(a);
2608}
2609
2610//--------------------------------------------------------------------------------------------------------------------------------------
2611
2612template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2613inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2J_dab(const index_t a, const index_t b) const
2614{
2615 return (a==b) ? 0.0 : 1.0 / ( m_data.mine().m_stretches(a) * m_data.mine().m_stretches(b) );
2616}
2617
2618//--------------------------------------------------------------------------------------------------------------------------------------
2619
2620template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2621inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_p() const
2622{
2623 return m_data.mine().m_stretches(2) * _dPsi_da(2);
2624}
2625
2626//--------------------------------------------------------------------------------------------------------------------------------------
2627
2628template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2629inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dp_da(const index_t a) const
2630{
2631 if (a==2)
2632 return m_data.mine().m_stretches(2) * _d2Psi_dab(2,a) + _dPsi_da(2);
2633 else
2634 return m_data.mine().m_stretches(2) * _d2Psi_dab(2,a);
2635
2636 // return m_data.mine().m_stretches(2) * _d2Psi_dab(2,a) + delta(a,2) * _dPsi_da(2);
2637}
2638
2639//--------------------------------------------------------------------------------------------------------------------------------------
2640
2641template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2642inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sa(const index_t a) const
2643{
2644 return _Sa_impl<comp>(a);
2645}
2646
2647template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2648template <bool _comp>
2649inline typename std::enable_if<_comp, T>::type
2650gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sa_impl(const index_t a) const
2651{
2652 return 1.0/m_data.mine().m_stretches(a) * _dPsi_da(a);
2653}
2654
2655template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2656template <bool _comp>
2657inline typename std::enable_if<!_comp, T>::type
2658gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sa_impl(const index_t a) const
2659{
2660 return 1.0/m_data.mine().m_stretches(a) * (_dPsi_da(a) - _p() * _dJ_da(a) );
2661}
2662
2663//--------------------------------------------------------------------------------------------------------------------------------------
2664
2665template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2666inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dSa_db(const index_t a, const index_t b) const
2667{
2668 return _dSa_db_impl<comp>(a,b);
2669}
2670
2671template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2672template <bool _comp>
2673inline typename std::enable_if<_comp, T>::type
2674gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dSa_db_impl(const index_t a, const index_t b) const
2675{
2676 T tmp = 1.0/m_data.mine().m_stretches(a) * _d2Psi_dab(a,b);
2677 if (a==b)
2678 tmp += - 1.0 / math::pow(m_data.mine().m_stretches(a),2) * _dPsi_da(a);
2679 return tmp;
2680}
2681
2682template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2683template <bool _comp>
2684inline typename std::enable_if<!_comp, T>::type
2685gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dSa_db_impl(const index_t a, const index_t b) const
2686{
2687 T tmp = 1.0/m_data.mine().m_stretches(a) * ( _d2Psi_dab(a,b) - _dp_da(a)*_dJ_da(b) - _dp_da(b)*_dJ_da(a) - _p() * _d2J_dab(a,b) );
2688 if (a==b)
2689 tmp += - 1.0 / math::pow(m_data.mine().m_stretches(a),2) * (_dPsi_da(a) - _p() * _dJ_da(a));
2690 return tmp;
2691}
2692
2693//--------------------------------------------------------------------------------------------------------------------------------------
2694
2695template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2696inline T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Cabcd(const index_t a, const index_t b, const index_t c, const index_t d) const
2697{
2698 return _Cabcd_impl<comp>(a,b,c,d);
2699}
2700
2701template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2702template <bool _comp>
2703inline typename std::enable_if<_comp, T>::type
2704gsMaterialMatrixNonlinear<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
2705{
2706 // Compute part with stress tensor involved.
2707 T frac = 0.0;
2708 T tmp = 0.0;
2709
2710 if (abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a)) < 1e-14)
2711 {
2712 // gsDebug<<"Stretches are equal; (abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a)) = "<<abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a))<<"\n";
2713 frac = 1.0 / (2.0 * m_data.mine().m_stretches(a) ) * ( _dSa_db(b,b) - _dSa_db(a,b));
2714 }
2715 else
2716 frac = ( _Sa(b)-_Sa(a) ) / (math::pow(m_data.mine().m_stretches(b),2) - math::pow(m_data.mine().m_stretches(a),2));
2717
2718 GISMO_ENSURE( ( (a < 3) && (b < 3) && (c < 3) && (d < 3) ) , "Index out of range. a="<<a<<", b="<<b<<", c="<<c<<", d="<<d);
2719 if ( ( (a==b) && (c==d)) )
2720 tmp = 1/m_data.mine().m_stretches(c) * _dSa_db(a,c);
2721 else if (( (a==d) && (b==c) && (a!=b) ) || ( ( (a==c) && (b==d) && (a!=b)) ))
2722 tmp = frac;
2723 // return 1/m_data.mine().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));
2724
2725 return tmp;
2726}
2727
2728template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2729template <bool _comp>
2730inline typename std::enable_if<!_comp, T>::type
2731gsMaterialMatrixNonlinear<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
2732{
2733 // Compute part with stress tensor involved.
2734 T frac = 0.0;
2735 T tmp = 0.0;
2736
2737 if (abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a)) < 1e-14)
2738 {
2739 // gsDebug<<"Stretches are equal; (abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a)) = "<<abs((m_data.mine().m_stretches(a) - m_data.mine().m_stretches(b)) / m_data.mine().m_stretches(a))<<"\n";
2740 frac = 1.0 / (2.0 * m_data.mine().m_stretches(a) ) * ( _dSa_db(b,b) - _dSa_db(a,b));
2741 }
2742 else
2743 frac = ( _Sa(b)-_Sa(a) ) / (math::pow(m_data.mine().m_stretches(b),2) - math::pow(m_data.mine().m_stretches(a),2));
2744
2745 GISMO_ENSURE( ( (a < 2) && (b < 2) && (c < 2) && (d < 2) ) , "Index out of range. a="<<a<<", b="<<b<<", c="<<c<<", d="<<d);
2746 if ( ( (a==b) && (c==d)) )
2747 tmp = 1/m_data.mine().m_stretches(c) * _dSa_db(a,c) + 1/(math::pow(m_data.mine().m_stretches(a),2) * math::pow(m_data.mine().m_stretches(c),2)) * ( math::pow(m_data.mine().m_stretches(2),2) * _d2Psi_dab(2,2) + 2*_dPsi_da(2)*m_data.mine().m_stretches(2) );
2748 else if (( (a==d) && (b==c) && (a!=b) ) || ( ( (a==c) && (b==d) && (a!=b)) ))
2749 tmp = frac;
2750 // return 1/m_data.mine().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))
2751 // + delta(a,b)*delta(c,d)*1/(math::pow(m_data.mine().m_stretches(a),2) * math::pow(m_data.mine().m_stretches(c),2)) * ( math::pow(m_data.mine().m_stretches(2),2) * _d2Psi_dab(2,2) + 2*_dPsi_da(2)*m_data.mine().m_stretches(2) );
2752
2753 return tmp;
2754}
2755
2756//--------------------------------------------------------------------------------------------------------------------------------------
2757
2758// template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2759// T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_dPsi(const index_t a) const
2760// {
2761// T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2762
2763// GISMO_ENSURE( ( (a < 3) ) , "Index out of range. a="<<a);
2764// GISMO_ENSURE(!comp,"Material model is not incompressible?");
2765
2766// if (m_material==9)
2767// {
2768// return mu * m_data.mine().m_stretches(a,0);
2769// }
2770// else if (m_material==3)
2771// {
2772// // return 2.0*m_data.mine().m_parvals.at(0)*m_data.mine().m_J0*m_data.mine().m_G(i,j)*m_data.mine().m_G(k,l) + m_data.mine().m_G(i,k)*m_data.mine().m_G(j,l) + m_data.mine().m_G(i,l)*m_data.mine().m_G(j,k);
2773// }
2774// else if (m_material==4)
2775// {
2776// gsMatrix<T> C(3,3);
2777// C.setZero();
2778// C.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
2779// C(2,2) = math::pow(m_J0,-2.0);
2780// this->_computeStretch(C,m_data.mine().m_gcon_ori);
2781
2782// return mu*m_data.mine().m_stretches.at(a);
2783// }
2784// else
2785// GISMO_ERROR("Material model not implemented.");
2786// }
2787
2788// template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2789// T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_d2Psi(const index_t a, const index_t b) const
2790// {
2791// T mu = m_data.mine().m_parvals.at(0) / (2. * (1. + m_data.mine().m_parvals.at(1)));
2792
2793// GISMO_ENSURE( ( (a < 3) && (b < 3) ) , "Index out of range. a="<<a<<", b="<<b);
2794// GISMO_ENSURE(!comp,"Material model is not incompressible?");
2795
2796// if (m_material==9)
2797// {
2798// return ( a==b ? mu : 0.0 );
2799// }
2800// else if (m_material==3)
2801// {
2802// // return 2.0*m_data.mine().m_parvals.at(0)*m_data.mine().m_J0*m_data.mine().m_G(i,j)*m_data.mine().m_G(k,l) + m_data.mine().m_G(i,k)*m_data.mine().m_G(j,l) + m_data.mine().m_G(i,l)*m_data.mine().m_G(j,k);
2803// }
2804// else if (m_material==4)
2805// {
2806// if (a==b)
2807// {
2808// return mu;
2809// }
2810// else
2811// return 0.0;
2812// }
2813// else
2814// GISMO_ERROR("Material model not implemented.");
2815// }
2816
2817namespace internal
2818{
2819
2823template<short_t d, class T, bool comp>
2824class gsXml< gsMaterialMatrixNonlinear<d,T,11,comp> >
2825{
2826private:
2827 gsXml() { }
2828 typedef gsMaterialMatrixNonlinear<d,T,11,comp> Object;
2829
2830public:
2831 GSXML_COMMON_FUNCTIONS(gsMaterialMatrixNonlinear<TMPLA4(d,T,11,comp)>);
2832 static std::string tag () { return "MaterialMatrix"; }
2833 static std::string type ()
2834 {
2835 std::string comp_str = ((comp) ? "Compressible" : "Incompressible");
2836 return comp_str + "NH" + to_string(d);
2837 }
2838
2839 GSXML_GET_POINTER(Object);
2840
2841 static void get_into(gsXmlNode * node,Object & obj)
2842 {
2843 obj = getMaterialMatrixFromXml< Object >( node );
2844 }
2845
2846 static gsXmlNode * put (const Object & obj,
2847 gsXmlTree & data)
2848 {
2849 return putMaterialMatrixToXml< Object >( obj,data );
2850 }
2851};
2852
2856template<short_t d, class T, bool comp>
2857class gsXml< gsMaterialMatrixNonlinear<d,T,12,comp> >
2858{
2859private:
2860 gsXml() { }
2861 typedef gsMaterialMatrixNonlinear<d,T,12,comp> Object;
2862
2863public:
2864 GSXML_COMMON_FUNCTIONS(gsMaterialMatrixNonlinear<TMPLA4(d,T,12,comp)>);
2865 static std::string tag () { return "MaterialMatrix"; }
2866 static std::string type ()
2867 {
2868 std::string comp_str = ((comp) ? "Compressible" : "Incompressible");
2869 return comp_str + "NHe" + to_string(d);
2870 }
2871
2872 GSXML_GET_POINTER(Object);
2873
2874 static void get_into(gsXmlNode * node,Object & obj)
2875 {
2876 obj = getMaterialMatrixFromXml< Object >( node );
2877 }
2878
2879 static gsXmlNode * put (const Object & obj,
2880 gsXmlTree & data)
2881 {
2882 return putMaterialMatrixToXml< Object >( obj,data );
2883 }
2884};
2885
2889template<short_t d, class T, bool comp>
2890class gsXml< gsMaterialMatrixNonlinear<d,T,13,comp> >
2891{
2892private:
2893 gsXml() { }
2894 typedef gsMaterialMatrixNonlinear<d,T,13,comp> Object;
2895
2896public:
2897 GSXML_COMMON_FUNCTIONS(gsMaterialMatrixNonlinear<TMPLA4(d,T,13,comp)>);
2898 static std::string tag () { return "MaterialMatrix"; }
2899 static std::string type ()
2900 {
2901 std::string comp_str = ((comp) ? "Compressible" : "Incompressible");
2902 return comp_str + "MR" + to_string(d);
2903 }
2904
2905 GSXML_GET_POINTER(Object);
2906
2907 static void get_into(gsXmlNode * node,Object & obj)
2908 {
2909 obj = getMaterialMatrixFromXml< Object >( node );
2910 }
2911
2912 static gsXmlNode * put (const Object & obj,
2913 gsXmlTree & data)
2914 {
2915 return putMaterialMatrixToXml< Object >( obj,data );
2916 }
2917};
2918
2922template<short_t d, class T, bool comp>
2923class gsXml< gsMaterialMatrixNonlinear<d,T,34,comp> >
2924{
2925private:
2926 gsXml() { }
2927 typedef gsMaterialMatrixNonlinear<d,T,34,comp> Object;
2928
2929public:
2930 GSXML_COMMON_FUNCTIONS(gsMaterialMatrixNonlinear<TMPLA4(d,T,34,comp)>);
2931 static std::string tag () { return "MaterialMatrix"; }
2932 static std::string type ()
2933 {
2934 std::string comp_str = ((comp) ? "Compressible" : "Incompressible");
2935 return comp_str + "OG" + to_string(d);
2936 }
2937
2938 GSXML_GET_POINTER(Object);
2939
2940 static void get_into(gsXmlNode * node,Object & obj)
2941 {
2942 obj = getMaterialMatrixFromXml< Object >( node );
2943 }
2944
2945 static gsXmlNode * put (const Object & obj,
2946 gsXmlTree & data)
2947 {
2948 return putMaterialMatrixToXml< Object >( obj,data );
2949 }
2950};
2951
2952}// namespace internal
2953
2954
2955
2956} // end namespace
gismo::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition gsAsMatrix.h:32
gismo::gsFunctionSet
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition gsFunctionSet.h:219
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition gsFunction.h:60
gismo::gsFunction::domainDim
virtual short_t domainDim() const=0
Dimension of the (source) domain.
gismo::gsFunction::targetDim
virtual short_t targetDim() const
Dimension of the target space.
Definition gsFunctionSet.h:595
gismo::gsFunction::eval_into
virtual void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const =0
Evaluate the function at points u into result.
gismo::gsMaterialMatrixBaseDim
This class defines the base class for material matrices.
Definition gsMaterialMatrixBaseDim.h:37
gismo::gsMaterialMatrixBase
This class defines the base class for material matrices.
Definition gsMaterialMatrixBase.h:33
gismo::gsMaterialMatrixBase::setParameters
virtual void setParameters(const std::vector< function_ptr > &pars)
Sets the material parameters.
Definition gsMaterialMatrixBase.h:684
gismo::gsMaterialMatrixNonlinear
This class defines hyperelastic material matrices.
Definition gsMaterialMatrixNonlinear.h:47
gismo::gsMaterialMatrixNonlinear::_Sij_impl
std::enable_if< _mat==Material::SvK &&_imp==Implementation::Analytical, T >::type _Sij_impl(const index_t i, const index_t j) const
Specialization for incompressible Sij(i,j) for SvK materials implemented analytically.
Definition gsMaterialMatrixNonlinear.hpp:1627
gismo::gsMaterialMatrixNonlinear::eval3D_vector
gsMatrix< T > eval3D_vector(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:549
gismo::gsMaterialMatrixNonlinear::eval3D_matrix_C
gsMatrix< T > eval3D_matrix_C(const gsMatrix< T > &Cmat, const index_t patch, const gsVector< T > &u, const T z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:304
gismo::gsMaterialMatrixNonlinear::_dJ_da
T _dJ_da(const index_t a) const
First derivative of the compressibilty function w.r.t. the stretche .
Definition gsMaterialMatrixNonlinear.hpp:2605
gismo::gsMaterialMatrixNonlinear::_d2Psi_dab
T _d2Psi_dab(const index_t a, const index_t b) const
Second derivative of a strain energy density function w.r.t. the stretches and .
Definition gsMaterialMatrixNonlinear.hpp:2447
gismo::gsMaterialMatrixNonlinear::eval3D_CauchyStress
gsMatrix< T > eval3D_CauchyStress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:961
gismo::gsMaterialMatrixNonlinear::_dSa_db_impl
std::enable_if< _com, T >::type _dSa_db_impl(const index_t a, const index_t b) const
Specialization of _dSa_db(a,b) for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:2674
gismo::gsMaterialMatrixNonlinear::_dPsi_da_vol
T _dPsi_da_vol(const index_t a) const
First derivative of the volumetric part of a strain energy density function w.r.t....
Definition gsMaterialMatrixNonlinear.hpp:2436
gismo::gsMaterialMatrixNonlinear::_eval3D_Compressible_detF
gsMatrix< T > _eval3D_Compressible_detF(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evaluates the jacobian determinant for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:926
gismo::gsMaterialMatrixNonlinear::pstretch_into
void pstretch_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:209
gismo::gsMaterialMatrixNonlinear::_eval3D_Incompressible_matrix
gsMatrix< T > _eval3D_Incompressible_matrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evalluates the incompressible material matrix.
Definition gsMaterialMatrixNonlinear.hpp:1190
gismo::gsMaterialMatrixNonlinear::_Cabcd_impl
std::enable_if< _com, T >::type _Cabcd_impl(const index_t a, const index_t b, const index_t c, const index_t d) const
Specialization of _Cabcd(a,b,c,d) for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:2704
gismo::gsMaterialMatrixNonlinear::_p
T _p() const
Lagrange multiplier for incompressible materials.
Definition gsMaterialMatrixNonlinear.hpp:2621
gismo::gsMaterialMatrixNonlinear::_dPsi
T _dPsi(const index_t i, const index_t j) const
Provides the derivative of the incompressible strain energy density function w.r.t....
Definition gsMaterialMatrixNonlinear.hpp:2067
gismo::gsMaterialMatrixNonlinear::_d2Psi_impl
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
Implementation of _d2Psi(i,j) for NH materials.
Definition gsMaterialMatrixNonlinear.hpp:2120
gismo::gsMaterialMatrixNonlinear::_initialize
void _initialize()
Initializes the object.
Definition gsMaterialMatrixNonlinear.hpp:202
gismo::gsMaterialMatrixNonlinear::setYoungsModulus
void setYoungsModulus(const gsFunctionSet< T > &YoungsModulus) override
Sets the YoungsModulus.
Definition gsMaterialMatrixNonlinear.hpp:1032
gismo::gsMaterialMatrixNonlinear::_Cijkl
T _Cijkl(const index_t i, const index_t j, const index_t k, const index_t l) const
Returns an entry of the material tensor C for incompressible materials.
Definition gsMaterialMatrixNonlinear.hpp:1231
gismo::gsMaterialMatrixNonlinear::_dI_1
T _dI_1(const index_t i, const index_t j) const
Provides the derivative of the first invariant compressible materials w.r.t. component C_{ij} of the ...
Definition gsMaterialMatrixNonlinear.hpp:2155
gismo::gsMaterialMatrixNonlinear::_eval3D_Compressible_matrix
gsMatrix< T > _eval3D_Compressible_matrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evalluates the compressible material matrix.
Definition gsMaterialMatrixNonlinear.hpp:2012
gismo::gsMaterialMatrixNonlinear::_pstretch_into_impl
std::enable_if< _com, void >::type _pstretch_into_impl(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const
Implementation of pstretch_into, specialization for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:218
gismo::gsMaterialMatrixNonlinear::_eval3D_matrix_impl
std::enable_if< _mat==Material::SvK, gsMatrix< T > >::type _eval3D_matrix_impl(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Implementation of the 3D (in-plane+thickness) evaluator of the material matrix, specialization for Sv...
Definition gsMaterialMatrixNonlinear.hpp:521
gismo::gsMaterialMatrixNonlinear::_d2Psi
T _d2Psi(const index_t i, const index_t j, const index_t k, const index_t l) const
Provides the second (mixed) derivative of the incompressible strain energy density function w....
Definition gsMaterialMatrixNonlinear.hpp:2109
gismo::gsMaterialMatrixNonlinear::_dPsi_da
T _dPsi_da(const index_t a) const
First derivative of a strain energy density function w.r.t. the stretch .
Definition gsMaterialMatrixNonlinear.hpp:2321
gismo::gsMaterialMatrixNonlinear::eval3D_pstressDir
gsMatrix< T > eval3D_pstressDir(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:597
gismo::gsMaterialMatrixNonlinear::_eval3D_stress_impl
std::enable_if< _mat==Material::SvK, gsMatrix< T > >::type _eval3D_stress_impl(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Implementation of the 3D (in-plane+thickness) evaluator of the stress vector, specialization for SvK ...
Definition gsMaterialMatrixNonlinear.hpp:668
gismo::gsMaterialMatrixNonlinear::_Sij
T _Sij(const index_t i, const index_t j) const
Returns an entry of the stress tensor S for incompressible materials.
Definition gsMaterialMatrixNonlinear.hpp:1619
gismo::gsMaterialMatrixNonlinear::_pstretchDir_into_impl
std::enable_if< _com, void >::type _pstretchDir_into_impl(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const
Implementation of stretchDir_into, specialization for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:266
gismo::gsMaterialMatrixNonlinear::defaultOptions
void defaultOptions() override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:195
gismo::gsMaterialMatrixNonlinear::_Sa
T _Sa(const index_t a) const
Component of the stress.
Definition gsMaterialMatrixNonlinear.hpp:2642
gismo::gsMaterialMatrixNonlinear::eval3D_detF
gsMatrix< T > eval3D_detF(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:892
gismo::gsMaterialMatrixNonlinear::print
std::ostream & print(std::ostream &os) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:153
gismo::gsMaterialMatrixNonlinear::_Sa_impl
std::enable_if< _com, T >::type _Sa_impl(const index_t a) const
Specialization of _Sa(a) for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:2650
gismo::gsMaterialMatrixNonlinear::_Cabcd
T _Cabcd(const index_t a, const index_t b, const index_t c, const index_t d) const
The material matrix for stretch-based implementations.
Definition gsMaterialMatrixNonlinear.hpp:2696
gismo::gsMaterialMatrixNonlinear::_dPsi_vol
T _dPsi_vol(const index_t i, const index_t j, const gsMatrix< T > &c, const gsMatrix< T > &cinv) const
Provides the derivative of the volumetric part of the compressible strain energy density function w....
Definition gsMaterialMatrixNonlinear.hpp:2229
gismo::gsMaterialMatrixNonlinear::_dI_2
T _dI_2(const index_t i, const index_t j, const gsMatrix< T > &c, const gsMatrix< T > &cinv) const
Provides the derivative of the second invariant for compressible materials w.r.t. component C_{ij} of...
Definition gsMaterialMatrixNonlinear.hpp:2163
gismo::gsMaterialMatrixNonlinear::_Cijkl3D_impl
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) const
Specialization for compressible Cijkl3D(i,j,k,l,c,cinv) for SvK materials implemented analytically.
Definition gsMaterialMatrixNonlinear.hpp:1423
gismo::gsMaterialMatrixNonlinear::eval3D_CauchyVector
gsMatrix< T > eval3D_CauchyVector(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:561
gismo::gsMaterialMatrixNonlinear::_eval3D_Incompressible_stress
gsMatrix< T > _eval3D_Incompressible_stress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evalluates the incompressible stress vector.
Definition gsMaterialMatrixNonlinear.hpp:825
gismo::gsMaterialMatrixNonlinear::pstretchDir_into
void pstretchDir_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:257
gismo::gsMaterialMatrixNonlinear::_d2Psi_vol
T _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
Provides the second (mixed) derivative of the volumetric part of the compressible strain energy densi...
Definition gsMaterialMatrixNonlinear.hpp:2307
gismo::gsMaterialMatrixNonlinear::_eval3D_Incompressible_detF
gsMatrix< T > _eval3D_Incompressible_detF(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evaluates the jacobian determinant for compressible materials.
Definition gsMaterialMatrixNonlinear.hpp:953
gismo::gsMaterialMatrixNonlinear::_eval3D_detF_impl
std::enable_if< _mat==Material::SvK, gsMatrix< T > >::type _eval3D_detF_impl(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Implementation of jacobina determinant (detF)
Definition gsMaterialMatrixNonlinear.hpp:901
gismo::gsMaterialMatrixNonlinear::_Cijkl3D
T _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
Returns an entry of the material tensor C for compressible materials without static condensation.
Definition gsMaterialMatrixNonlinear.hpp:1414
gismo::gsMaterialMatrixNonlinear::_Cijkl_impl
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) const
Specialization for incompressible Cijkl(i,j,k,l) for SvK materials implemented analytically.
Definition gsMaterialMatrixNonlinear.hpp:1242
gismo::gsMaterialMatrixNonlinear::_dSa_db
T _dSa_db(const index_t a, const index_t b) const
First derivative of the component of the stress w.r.t. the stretch .
Definition gsMaterialMatrixNonlinear.hpp:2666
gismo::gsMaterialMatrixNonlinear::_d2Psi_dab_vol
T _d2Psi_dab_vol(const index_t a, const index_t b) const
Second derivative of the volumetric part of a strain energy density function w.r....
Definition gsMaterialMatrixNonlinear.hpp:2592
gismo::gsMaterialMatrixNonlinear::setPoissonsRatio
void setPoissonsRatio(const gsFunctionSet< T > &PoissonsRatio) override
Sets the Poisson's Ratio.
Definition gsMaterialMatrixNonlinear.hpp:1044
gismo::gsMaterialMatrixNonlinear::eval3D_stress_C
gsMatrix< T > eval3D_stress_C(const gsMatrix< T > &Cmat, const index_t patch, const gsVector< T > &u, const T z, enum MaterialOutput out=MaterialOutput::Generic) const
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:693
gismo::gsMaterialMatrixNonlinear::eval3D_stress
gsMatrix< T > eval3D_stress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:659
gismo::gsMaterialMatrixNonlinear::eval3D_vector_C
gsMatrix< T > eval3D_vector_C(const gsMatrix< T > &Cmat, const index_t patch, const gsVector< T > &u, const T z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:555
gismo::gsMaterialMatrixNonlinear::_eval3D_Compressible_C33
gsMatrix< T > _eval3D_Compressible_C33(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evaluates the C33 (through thickness) component of the deformation tensor C for compressible material...
Definition gsMaterialMatrixNonlinear.hpp:1815
gismo::gsMaterialMatrixNonlinear::getYoungsModulus
const function_ptr getYoungsModulus() const override
Gets the YoungsModulus.
Definition gsMaterialMatrixNonlinear.hpp:1038
gismo::gsMaterialMatrixNonlinear::_Sii
T _Sii(const index_t i) const
Returns an entry of the diagonal of the stress tensor S for incompressible materials.
Definition gsMaterialMatrixNonlinear.hpp:1798
gismo::gsMaterialMatrixNonlinear::gsMaterialMatrixNonlinear
gsMaterialMatrixNonlinear()
Destructor.
Definition gsMaterialMatrixNonlinear.h:124
gismo::gsMaterialMatrixNonlinear::getPoissonsRatio
const function_ptr getPoissonsRatio() const override
Gets the Poisson's Ratio.
Definition gsMaterialMatrixNonlinear.hpp:1050
gismo::gsMaterialMatrixNonlinear::_dp_da
T _dp_da(const index_t a) const
First derivative of the Lagrange multiplier for incompressible materials w.r.t. the stretch .
Definition gsMaterialMatrixNonlinear.hpp:2629
gismo::gsMaterialMatrixNonlinear::_d2J_dab
T _d2J_dab(const index_t a, const index_t b) const
First derivative of the compressibilty function w.r.t. the stretches and .
Definition gsMaterialMatrixNonlinear.hpp:2613
gismo::gsMaterialMatrixNonlinear::eval3D_pstress
gsMatrix< T > eval3D_pstress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:567
gismo::gsMaterialMatrixNonlinear::eval3D_dmatrix
gsMatrix< T > eval3D_dmatrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:344
gismo::gsMaterialMatrixNonlinear::_eval3D_Compressible_stress
gsMatrix< T > _eval3D_Compressible_stress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const
Evalluates the compressible stress vector.
Definition gsMaterialMatrixNonlinear.hpp:727
gismo::gsMaterialMatrixNonlinear::_dPsi_impl
std::enable_if< _mat==Material::NH, T >::type _dPsi_impl(const index_t i, const index_t j) const
Implementation of _dPsi(i,j) for NH materials.
Definition gsMaterialMatrixNonlinear.hpp:2078
gismo::gsMaterialMatrixNonlinear::eval3D_matrix
gsMatrix< T > eval3D_matrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:509
gismo::gsMaterialMatrixNonlinear::eval3D_CauchyPStress
gsMatrix< T > eval3D_CauchyPStress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixNonlinear.hpp:627
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshapeCol
gsAsMatrix< T, Dynamic, Dynamic > reshapeCol(index_t c, index_t n, index_t m)
Returns column c of the matrix resized to n x m matrix This function assumes that the matrix is size ...
Definition gsMatrix.h:231
gismo::gsMatrix::reshape
gsAsMatrix< T, Dynamic, Dynamic > reshape(index_t n, index_t m)
Returns the matrix resized to n x m matrix (data is not copied) This function assumes that the matrix...
Definition gsMatrix.h:221
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition gsMaterialMatrixUtils.h:99
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsFunction.h
Provides declaration of Function abstract interface.
gsMaterialMatrixNonlinear.h
Provides hyperelastic material matrices.
example_shell3D.c1
c1
See gismo/filedata/surfaces/simple.xml for the geometry.
Definition example_shell3D.py:40
gismo::internal::to_string
std::string to_string(const unsigned &i)
Helper to convert small unsigned to string.
Definition gsXml.cpp:74
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::gsAllCloseRelativeToMax
bool gsAllCloseRelativeToMax(const matrix_t1 &a, const matrix_t2 &b, const typename matrix_t1::Scalar &tol)
tests if the difference between two matrices is bounded by tol in norm
Definition gsMath.h:452