G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
28 using namespace gismo;
29 
30 template <short_t d, typename T>
31 class Cfun : public gsFunction<T>
32 {
33 public:
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 
76 private:
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 
83 namespace gismo
84 {
85 
86 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
87 gsMaterialMatrixNonlinear<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 
97 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
98 gsMaterialMatrixNonlinear<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 
107 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
108 gsMaterialMatrixNonlinear<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 
118 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
119 gsMaterialMatrixNonlinear<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 
127 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
128 gsMaterialMatrixNonlinear<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 
137 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
138 gsMaterialMatrixNonlinear<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 
152 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
153 std::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 
194 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
195 void 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 
201 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
202 void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_initialize()
203 {
204  // Set default options
205  this->defaultOptions();
206 }
207 
208 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
209 void 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 
215 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
216 template <bool _comp>
217 typename std::enable_if<!_comp, void>::type
218 gsMaterialMatrixNonlinear<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 
223 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
224 template <bool _comp>
225 typename std::enable_if<_comp, void>::type
226 gsMaterialMatrixNonlinear<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 
256 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
257 void 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 
263 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
264 template <bool _comp>
265 typename std::enable_if<!_comp, void>::type
266 gsMaterialMatrixNonlinear<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 
271 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
272 template <bool _comp>
273 typename std::enable_if<_comp, void>::type
274 gsMaterialMatrixNonlinear<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 
303 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
304 gsMatrix<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 
313 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
314 template <enum Material _mat, bool _comp>
315 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
316 gsMaterialMatrixNonlinear<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 
323 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
324 template <enum Material _mat, bool _comp>
325 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
326 gsMaterialMatrixNonlinear<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 
333 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
334 template <enum Material _mat, bool _comp>
335 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
336 gsMaterialMatrixNonlinear<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 
343 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
344 gsMatrix<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 
371 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
372 constexpr 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 
377 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
378 template <enum Material _mat, bool _comp>
379 constexpr typename std::enable_if<!(!_comp && _mat==Material::NH), gsMatrix<T>>::type
380 gsMaterialMatrixNonlinear<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 
394 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
395 template <enum Material _mat, bool _comp>
396 constexpr typename std::enable_if<(!_comp && _mat==Material::NH), gsMatrix<T>>::type
397 gsMaterialMatrixNonlinear<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 
438 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
439 constexpr 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 
451 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
452 constexpr 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 
457 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
458 template <enum Material _mat, bool _comp>
459 constexpr typename std::enable_if<(!_comp && _mat==Material::NH), T>::type
460 gsMaterialMatrixNonlinear<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 
500 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
501 template <enum Material _mat, bool _comp>
502 constexpr typename std::enable_if<!(!_comp && _mat==Material::NH), T>::type
503 gsMaterialMatrixNonlinear<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 
508 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
509 gsMatrix<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 
518 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
519 template <enum Material _mat, bool _comp>
520 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
521 gsMaterialMatrixNonlinear<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 
528 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
529 template <enum Material _mat, bool _comp>
530 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
531 gsMaterialMatrixNonlinear<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 
538 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
539 template <enum Material _mat, bool _comp>
540 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
541 gsMaterialMatrixNonlinear<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 
548 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
549 gsMatrix<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 
554 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
555 gsMatrix<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 
560 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
561 gsMatrix<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 
566 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
567 gsMatrix<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 
596 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
597 gsMatrix<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 
626 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
627 gsMatrix<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 
658 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
659 gsMatrix<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 
665 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
666 template <enum Material _mat, bool _comp>
667 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
668 gsMaterialMatrixNonlinear<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 
674 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
675 template <enum Material _mat, bool _comp>
676 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
677 gsMaterialMatrixNonlinear<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 
683 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
684 template <enum Material _mat, bool _comp>
685 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
686 gsMaterialMatrixNonlinear<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 
692 template <short_t dim, class T, index_t matId, bool comp, enum Material mat, enum Implementation imp >
693 gsMatrix<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 
699 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
700 template <enum Material _mat, bool _comp>
701 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
702 gsMaterialMatrixNonlinear<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 
708 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
709 template <enum Material _mat, bool _comp>
710 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
711 gsMaterialMatrixNonlinear<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 
717 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
718 template <enum Material _mat, bool _comp>
719 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
720 gsMaterialMatrixNonlinear<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 
726 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
727 gsMatrix<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 
775 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
776 gsMatrix<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 
824 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
825 gsMatrix<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 
860 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
861 gsMatrix<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 
891 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
892 gsMatrix<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 
898 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
899 template <enum Material _mat, bool _comp>
900 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
901 gsMaterialMatrixNonlinear<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 
907 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
908 template <enum Material _mat, bool _comp>
909 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
910 gsMaterialMatrixNonlinear<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 
916 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
917 template <enum Material _mat, bool _comp>
918 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
919 gsMaterialMatrixNonlinear<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 
925 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
926 gsMatrix<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 
952 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
953 gsMatrix<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 
960 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
961 gsMatrix<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 
967 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
968 template <enum Material _mat, bool _comp>
969 typename std::enable_if<_mat==Material::SvK, gsMatrix<T>>::type
970 gsMaterialMatrixNonlinear<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 
976 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
977 template <enum Material _mat, bool _comp>
978 typename std::enable_if<!_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
979 gsMaterialMatrixNonlinear<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 
985 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
986 template <enum Material _mat, bool _comp>
987 typename std::enable_if<_comp && !(_mat==Material::SvK), gsMatrix<T>>::type
988 gsMaterialMatrixNonlinear<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 
994 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
995 gsMatrix<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 
1024 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1025 gsMatrix<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 
1031 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1032 void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::setYoungsModulus(const gsFunctionSet<T> & YoungsModulus)
1033 {
1034  Base::setParameter(0,YoungsModulus);
1035 }
1036 
1037 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1038 const 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 
1043 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1044 void gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::setPoissonsRatio(const gsFunctionSet<T> & PoissonsRatio)
1045 {
1046  Base::setParameter(1,PoissonsRatio);
1047 }
1048 
1049 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1050 const 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 
1055 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1056 template <enum Material _mat>
1057 typename std::enable_if<_mat==Material::MR, void>::type
1058 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setRatio_impl(const gsFunctionSet<T> & Ratio)
1059 {
1060  Base::setParameter(2,Ratio);
1061 }
1062 
1063 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1064 template <enum Material _mat>
1065 typename std::enable_if<_mat!=Material::MR, void>::type
1066 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setRatio_impl(const gsFunctionSet<T> & Ratio)
1067 {
1068  GISMO_NO_IMPLEMENTATION;
1069 }
1070 
1071 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1072 template <enum Material _mat>
1073 typename std::enable_if<_mat==Material::MR, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1074 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getRatio_impl() const
1075 {
1076  return Base::getParameter(2);
1077 }
1078 
1079 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1080 template <enum Material _mat>
1081 typename std::enable_if<_mat!=Material::MR, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1082 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getRatio_impl() const
1083 {
1084  GISMO_NO_IMPLEMENTATION;
1085 }
1086 
1087 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1088 template <enum Material _mat>
1089 typename std::enable_if<_mat==Material::OG, void>::type
1090 gsMaterialMatrixNonlinear<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 
1095 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1096 template <enum Material _mat>
1097 typename std::enable_if<_mat!=Material::OG, void>::type
1098 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setMu_impl(const index_t & i, const gsFunctionSet<T> & Mu_i)
1099 {
1100  GISMO_NO_IMPLEMENTATION;
1101 }
1102 
1103 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1104 template <enum Material _mat>
1105 typename std::enable_if<_mat==Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1106 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getMu_impl(const index_t & i) const
1107 {
1108  return Base::getParameter(2+2*i);
1109 }
1110 
1111 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1112 template <enum Material _mat>
1113 typename std::enable_if<_mat!=Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1114 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getMu_impl(const index_t & i) const
1115 {
1116  GISMO_NO_IMPLEMENTATION;
1117 }
1118 
1119 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1120 template <enum Material _mat>
1121 typename std::enable_if<_mat==Material::OG, void>::type
1122 gsMaterialMatrixNonlinear<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 
1127 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1128 template <enum Material _mat>
1129 typename std::enable_if<_mat!=Material::OG, void>::type
1130 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_setAlpha_impl(const index_t & i, const gsFunctionSet<T> & Alpha_i)
1131 {
1132  GISMO_NO_IMPLEMENTATION;
1133 }
1134 
1135 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1136 template <enum Material _mat>
1137 typename std::enable_if<_mat==Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1138 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getAlpha_impl(const index_t & i) const
1139 {
1140  return Base::getParameter(2+2*i+1);
1141 }
1142 
1143 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1144 template <enum Material _mat>
1145 typename std::enable_if<_mat!=Material::OG, const typename gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::function_ptr >::type
1146 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_getAlpha_impl(const index_t & i) const
1147 {
1148  GISMO_NO_IMPLEMENTATION;
1149 }
1150 
1151 
1152 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1153 gsMatrix<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 
1189 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1190 gsMatrix<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 
1230 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1231 constexpr 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 
1239 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1240 template <enum Material _mat, enum Implementation _imp>
1241 constexpr typename std::enable_if<_mat==Material::SvK && _imp==Implementation::Analytical, T>::type
1242 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
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 
1257 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1258 template <enum Material _mat, enum Implementation _imp>
1259 constexpr typename std::enable_if<_mat==Material::NH && _imp==Implementation::Analytical, T>::type
1260 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
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 
1269 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1270 template <enum Material _mat, enum Implementation _imp>
1271 constexpr typename std::enable_if<_mat==Material::MR && _imp==Implementation::Analytical, T>::type
1272 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
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 
1295 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1296 template <enum Material _mat, enum Implementation _imp>
1297 constexpr typename std::enable_if<_imp==Implementation::Spectral, T>::type
1298 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
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 
1352 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1353 template <enum Material _mat, enum Implementation _imp>
1354 constexpr typename std::enable_if<_imp==Implementation::Generalized, T>::type
1355 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
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
1381 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1382 constexpr 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 
1394 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1395 template <enum Implementation _imp>
1396 constexpr typename std::enable_if< _imp==Implementation::Spectral, T>::type
1397 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 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 
1404 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1405 template <enum Implementation _imp>
1406 constexpr typename std::enable_if<!(_imp==Implementation::Spectral), T>::type
1407 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 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
1413 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1414 constexpr 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 
1420 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1421 template <enum Material _mat, enum Implementation _imp>
1422 constexpr typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
1423 gsMaterialMatrixNonlinear<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 
1428 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1429 template <enum Material _mat, enum Implementation _imp>
1430 constexpr typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1431 gsMaterialMatrixNonlinear<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 
1452 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1453 template <enum Material _mat, enum Implementation _imp>
1454 constexpr typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1455 gsMaterialMatrixNonlinear<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 
1487 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1488 template <enum Material _mat, enum Implementation _imp>
1489 constexpr typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
1490 gsMaterialMatrixNonlinear<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 
1503 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1504 template <enum Material _mat, enum Implementation _imp>
1505 constexpr typename std::enable_if<_imp == Implementation::Spectral, T>::type
1506 gsMaterialMatrixNonlinear<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 
1567 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1568 template <enum Material _mat, enum Implementation _imp>
1569 constexpr typename std::enable_if<_imp == Implementation::Generalized, T>::type
1570 gsMaterialMatrixNonlinear<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 
1618 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1619 constexpr 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 
1624 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1625 template <enum Material _mat, enum Implementation _imp>
1626 constexpr typename std::enable_if<_mat==Material::SvK && _imp == Implementation::Analytical, T>::type
1627 gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sij_impl(const index_t i, const index_t j) const
1628 {
1629  GISMO_ERROR("Not implemented");
1630 }
1631 
1632 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1633 template <enum Material _mat, enum Implementation _imp>
1634 constexpr typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1635 gsMaterialMatrixNonlinear<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 
1644 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1645 template <enum Material _mat, enum Implementation _imp>
1646 constexpr typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1647 gsMaterialMatrixNonlinear<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 
1667 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1668 template <enum Material _mat, enum Implementation _imp>
1669 constexpr typename std::enable_if<_imp == Implementation::Spectral, T>::type
1670 gsMaterialMatrixNonlinear<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 
1689 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1690 template <enum Material _mat, enum Implementation _imp>
1691 constexpr typename std::enable_if<_imp == Implementation::Generalized, T>::type
1692 gsMaterialMatrixNonlinear<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 
1702 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1703 constexpr 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 
1708 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1709 template <enum Material _mat, enum Implementation _imp>
1710 constexpr typename std::enable_if<_mat==Material::NH && _imp == Implementation::Analytical, T>::type
1711 gsMaterialMatrixNonlinear<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 
1728 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1729 template <enum Material _mat, enum Implementation _imp>
1730 constexpr typename std::enable_if<_mat==Material::MR && _imp == Implementation::Analytical, T>::type
1731 gsMaterialMatrixNonlinear<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 
1755 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1756 template <enum Material _mat, enum Implementation _imp>
1757 constexpr typename std::enable_if<_mat==Material::NH_ext && _imp == Implementation::Analytical, T>::type
1758 gsMaterialMatrixNonlinear<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 
1768 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1769 template <enum Material _mat, enum Implementation _imp>
1770 constexpr typename std::enable_if<_imp == Implementation::Spectral, T>::type
1771 gsMaterialMatrixNonlinear<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 
1784 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1785 template <enum Material _mat, enum Implementation _imp>
1786 constexpr typename std::enable_if<_imp == Implementation::Generalized, T>::type
1787 gsMaterialMatrixNonlinear<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 
1797 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1798 constexpr 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 
1807 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1808 constexpr 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 
1814 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1815 gsMatrix<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 
1888 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1889 gsMatrix<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 
1961 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
1962 gsMatrix<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 
2011 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2012 gsMatrix<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 // ---------------------------------------------------------------------------------------------------------------------------------
2066 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2067 constexpr 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 
2075 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2076 template <enum Material _mat>
2077 constexpr typename std::enable_if<_mat==Material::NH, T>::type
2078 gsMaterialMatrixNonlinear<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 
2084 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2085 template <enum Material _mat>
2086 constexpr typename std::enable_if<_mat==Material::MR, T>::type
2087 gsMaterialMatrixNonlinear<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 
2108 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2109 constexpr 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 
2117 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2118 template <enum Material _mat>
2119 constexpr typename std::enable_if<_mat==Material::NH, T>::type
2120 gsMaterialMatrixNonlinear<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 
2125 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2126 template <enum Material _mat>
2127 constexpr typename std::enable_if<_mat==Material::MR, T>::type
2128 gsMaterialMatrixNonlinear<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 
2154 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2155 constexpr 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 
2162 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2163 constexpr 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 
2174 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2175 constexpr 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 
2183 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2184 template <enum Material _mat>
2185 constexpr typename std::enable_if<_mat==Material::NH, T>::type
2186 gsMaterialMatrixNonlinear<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 
2197 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2198 template <enum Material _mat>
2199 constexpr typename std::enable_if<_mat==Material::MR, T>::type
2200 gsMaterialMatrixNonlinear<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 
2216 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2217 template <enum Material _mat>
2218 constexpr typename std::enable_if<_mat==Material::NH_ext, T>::type
2219 gsMaterialMatrixNonlinear<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.
2228 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2229 constexpr 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 
2237 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2238 constexpr 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 
2246 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2247 template <enum Material _mat>
2248 constexpr typename std::enable_if<_mat==Material::NH, T>::type
2249 gsMaterialMatrixNonlinear<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 
2265 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2266 template <enum Material _mat>
2267 constexpr typename std::enable_if<_mat==Material::MR, T>::type
2268 gsMaterialMatrixNonlinear<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 
2295 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2296 template <enum Material _mat>
2297 constexpr typename std::enable_if<_mat==Material::NH_ext, T>::type
2298 gsMaterialMatrixNonlinear<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 
2306 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2307 constexpr 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 
2320 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2321 constexpr 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 
2327 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2328 template <enum Material _mat, bool _comp>
2329 constexpr typename std::enable_if<_comp && (_mat==Material::NH), T>::type
2330 gsMaterialMatrixNonlinear<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 
2341 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2342 template <enum Material _mat, bool _comp>
2343 constexpr typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
2344 gsMaterialMatrixNonlinear<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 
2351 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2352 template <enum Material _mat, bool _comp>
2353 constexpr typename std::enable_if<_comp && (_mat==Material::MR), T>::type
2354 gsMaterialMatrixNonlinear<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 
2371 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2372 template <enum Material _mat, bool _comp>
2373 constexpr typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
2374 gsMaterialMatrixNonlinear<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 
2386 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2387 template <enum Material _mat, bool _comp>
2388 constexpr typename std::enable_if<_comp && (_mat==Material::OG), T>::type
2389 gsMaterialMatrixNonlinear<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 }
2404 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2405 template <enum Material _mat, bool _comp>
2406 constexpr typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
2407 gsMaterialMatrixNonlinear<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 
2421 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2422 template <enum Material _mat, bool _comp>
2423 constexpr typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
2424 gsMaterialMatrixNonlinear<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 
2435 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2436 constexpr 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 
2446 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2447 constexpr 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 
2453 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2454 template <enum Material _mat, bool _comp>
2455 constexpr typename std::enable_if<_comp && (_mat==Material::NH), T>::type
2456 gsMaterialMatrixNonlinear<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 
2473 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2474 template <enum Material _mat, bool _comp>
2475 constexpr typename std::enable_if<!_comp && (_mat==Material::NH), T>::type
2476 gsMaterialMatrixNonlinear<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 
2483 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2484 template <enum Material _mat, bool _comp>
2485 constexpr typename std::enable_if<_comp && (_mat==Material::MR), T>::type
2486 gsMaterialMatrixNonlinear<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 
2520 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2521 template <enum Material _mat, bool _comp>
2522 constexpr typename std::enable_if<!_comp && (_mat==Material::MR), T>::type
2523 gsMaterialMatrixNonlinear<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 
2539 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2540 template <enum Material _mat, bool _comp>
2541 constexpr typename std::enable_if<_comp && (_mat==Material::OG), T>::type
2542 gsMaterialMatrixNonlinear<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 
2562 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2563 template <enum Material _mat, bool _comp>
2564 constexpr typename std::enable_if<!_comp && (_mat==Material::OG), T>::type
2565 gsMaterialMatrixNonlinear<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 
2579 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2580 template <enum Material _mat, bool _comp>
2581 constexpr typename std::enable_if<_comp && (_mat==Material::NH_ext), T>::type
2582 gsMaterialMatrixNonlinear<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 
2591 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2592 constexpr 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 
2604 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2605 constexpr 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 
2612 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2613 constexpr 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 
2620 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2621 constexpr 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 
2628 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2629 constexpr 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 
2641 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2642 constexpr T gsMaterialMatrixNonlinear<dim,T,matId,comp,mat,imp>::_Sa(const index_t a) const
2643 {
2644  return _Sa_impl<comp>(a);
2645 }
2646 
2647 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2648 template <bool _comp>
2649 constexpr typename std::enable_if<_comp, T>::type
2650 gsMaterialMatrixNonlinear<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 
2655 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2656 template <bool _comp>
2657 constexpr typename std::enable_if<!_comp, T>::type
2658 gsMaterialMatrixNonlinear<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 
2665 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2666 constexpr 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 
2671 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2672 template <bool _comp>
2673 constexpr typename std::enable_if<_comp, T>::type
2674 gsMaterialMatrixNonlinear<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 
2682 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2683 template <bool _comp>
2684 constexpr typename std::enable_if<!_comp, T>::type
2685 gsMaterialMatrixNonlinear<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 
2695 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2696 constexpr 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 
2701 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2702 template <bool _comp>
2703 constexpr typename std::enable_if<_comp, T>::type
2704 gsMaterialMatrixNonlinear<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 
2728 template <short_t dim, class T, short_t matId, bool comp, enum Material mat, enum Implementation imp >
2729 template <bool _comp>
2730 constexpr typename std::enable_if<!_comp, T>::type
2731 gsMaterialMatrixNonlinear<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 
2817 namespace internal
2818 {
2819 
2823 template<short_t d, class T, bool comp>
2824 class gsXml< gsMaterialMatrixNonlinear<d,T,11,comp> >
2825 {
2826 private:
2827  gsXml() { }
2828  typedef gsMaterialMatrixNonlinear<d,T,11,comp> Object;
2829 
2830 public:
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 
2856 template<short_t d, class T, bool comp>
2857 class gsXml< gsMaterialMatrixNonlinear<d,T,12,comp> >
2858 {
2859 private:
2860  gsXml() { }
2861  typedef gsMaterialMatrixNonlinear<d,T,12,comp> Object;
2862 
2863 public:
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 
2889 template<short_t d, class T, bool comp>
2890 class gsXml< gsMaterialMatrixNonlinear<d,T,13,comp> >
2891 {
2892 private:
2893  gsXml() { }
2894  typedef gsMaterialMatrixNonlinear<d,T,13,comp> Object;
2895 
2896 public:
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 
2922 template<short_t d, class T, bool comp>
2923 class gsXml< gsMaterialMatrixNonlinear<d,T,34,comp> >
2924 {
2925 private:
2926  gsXml() { }
2927  typedef gsMaterialMatrixNonlinear<d,T,34,comp> Object;
2928 
2929 public:
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::gsMaterialMatrixNonlinear::_d2Psi_vol
constexpr 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::MaterialOutput
MaterialOutput
This class describes the output type.
Definition: gsMaterialMatrixUtils.h:98
gismo::gsMaterialMatrixNonlinear::setPoissonsRatio
void setPoissonsRatio(const gsFunctionSet< T > &PoissonsRatio) override
Sets the Poisson&#39;s Ratio.
Definition: gsMaterialMatrixNonlinear.hpp:1044
gismo::gsMaterialMatrixNonlinear::_Cijkl3D_impl
constexpr 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::_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.
gismo::gsMaterialMatrixBase::setParameters
virtual void setParameters(const std::vector< function_ptr > &pars)
Sets the material parameters.
Definition: gsMaterialMatrixBase.h:684
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
This class defines hyperelastic material matrices.
Definition: gsMaterialMatrixNonlinear.h:47
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
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::_dPsi_impl
constexpr 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_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::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition: gsLinearAlgebra.h:126
gismo::gsMaterialMatrixNonlinear::_dp_da
constexpr 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
short_t
#define short_t
Definition: gsConfig.h:35
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::defaultOptions
void defaultOptions() override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixNonlinear.hpp:195
gismo::gsMaterialMatrixNonlinear::_dSa_db_impl
constexpr 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.
gismo::gsMaterialMatrixNonlinear::_Cijkl3D
constexpr 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::_dPsi_da
constexpr 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_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_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::_Cabcd
constexpr 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::_Sa
constexpr T _Sa(const index_t a) const
Component of the stress.
Definition: gsMaterialMatrixNonlinear.hpp:2642
gismo::gsMaterialMatrixNonlinear::_Cijkl
constexpr 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
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
gismo::gsMaterialMatrixNonlinear::_d2Psi_dab
constexpr 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::setYoungsModulus
void setYoungsModulus(const gsFunctionSet< T > &YoungsModulus) override
Sets the YoungsModulus.
Definition: gsMaterialMatrixNonlinear.hpp:1032
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition: gsFunction.h:59
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::_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::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_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gismo::gsMaterialMatrixBase
This class defines the base class for material matrices.
Definition: gsMaterialMatrixBase.h:32
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::gsMaterialMatrixNonlinear::_dI_2
constexpr 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::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition: gsMatrix.h:211
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::_Cijkl_impl
constexpr 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::_dPsi_vol
constexpr 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::internal::to_string
std::string to_string(const unsigned &i)
Helper to convert small unsigned to string.
Definition: gsXml.cpp:74
gsWarn
#define gsWarn
Definition: gsDebug.h:50
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::getPoissonsRatio
const function_ptr getPoissonsRatio() const override
Gets the Poisson&#39;s Ratio.
Definition: gsMaterialMatrixNonlinear.hpp:1050
gismo::gsMaterialMatrixNonlinear::_d2J_dab
constexpr 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::_dSa_db
constexpr 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::_Sij_impl
constexpr 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::_d2Psi
constexpr 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::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::_p
constexpr T _p() const
Lagrange multiplier for incompressible materials.
Definition: gsMaterialMatrixNonlinear.hpp:2621
gismo::gsMaterialMatrixNonlinear::_Cabcd_impl
constexpr 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.
gismo::gsMaterialMatrixNonlinear::_d2Psi_impl
constexpr 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::gsFunctionSet
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
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::gsMaterialMatrixNonlinear::getYoungsModulus
const function_ptr getYoungsModulus() const override
Gets the YoungsModulus.
Definition: gsMaterialMatrixNonlinear.hpp:1038
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::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::gsMaterialMatrixNonlinear::_Sii
constexpr 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:125
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
gsMaterialMatrixNonlinear.h
Provides hyperelastic material matrices.
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::_dI_1
constexpr 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::_initialize
void _initialize()
Initializes the object.
Definition: gsMaterialMatrixNonlinear.hpp:202
gismo::gsMaterialMatrixNonlinear::_dPsi_da_vol
constexpr T _dPsi_da_vol(const index_t a) const
First derivative of the volumetric part of a strain energy density function w.r.t. the stretch .
Definition: gsMaterialMatrixNonlinear.hpp:2436
gismo::gsMaterialMatrixNonlinear::_d2Psi_dab_vol
constexpr 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.t. the stretches and .
Definition: gsMaterialMatrixNonlinear.hpp:2592
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::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::_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::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_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
gismo::gsMaterialMatrixNonlinear::_d2Psi_dab_impl
constexpr std::enable_if< _com &&(_mat==Material::NH), T >::type _d2Psi_dab_impl(const index_t a, const index_t b) const
Specialization of _d2Psi_dab(a,b) for compressible NH materials.
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
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.
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_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::_dJ_da
constexpr 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::_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
gsFunction.h
Provides declaration of Function abstract interface.
gismo::gsMaterialMatrixNonlinear::print
std::ostream & print(std::ostream &os) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixNonlinear.hpp:153
gismo::gsMaterialMatrixNonlinear::_dPsi_da_impl
constexpr std::enable_if< _com &&(_mat==Material::NH), T >::type _dPsi_da_impl(const index_t a) const
Specialization of _dPsi_da(a) for compressible NH materials.
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::_Sij
constexpr 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::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
gismo::gsMaterialMatrixNonlinear::_dPsi
constexpr T _dPsi(const index_t i, const index_t j) const
Provides the derivative of the incompressible strain energy density function w.r.t. component C_{ij} of the deformation tensor.
Definition: gsMaterialMatrixNonlinear.hpp:2067
gismo::gsMaterialMatrixNonlinear::_Sa_impl
constexpr std::enable_if< _com, T >::type _Sa_impl(const index_t a) const
Specialization of _Sa(a) for compressible materials.