G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsMaterialMatrixBaseDim.hpp
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 
23 
24 #pragma once
25 
26 #include <gsKLShell/src/gsMaterialMatrixBaseDim.h>
27 #include <gsCore/gsGeometry.h>
28 
29 #include <gsKLShell/src/gsMaterialMatrixLinear.h>
30 #include <gsKLShell/src/gsMaterialMatrixNonlinear.h>
31 
32 namespace gismo
33 {
34 
35 //--------------------------------------------------------------------------------------------------------------------------------------
36 
37 template <short_t dim, class T >
38 void gsMaterialMatrixBaseDim<dim,T>::defaultOptions()
39 {
40  m_options.addInt("TensionField","Tension field definition - 0: Strain-based, 1: Stress-based, 2: Mixed",2);
41 }
42 
43 template <short_t dim, class T >
44 void gsMaterialMatrixBaseDim<dim,T>::density_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
45 {
46  gsMapData<T> map;
47  map.flags = NEED_VALUE;
48  map.points = u;
49  static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
50 
51  result.resize(1, u.cols());
52  m_thickness->piece(patch).eval_into(map.values[0], m_data.mine().m_Tmat);
53  m_density->piece(patch).eval_into(map.values[0], m_data.mine().m_rhomat);
54  for (index_t i = 0; i != u.cols(); ++i) // points
55  {
56  result(0,i) = m_data.mine().m_Tmat(0,i)*m_data.mine().m_rhomat(0,i);
57  }
58 }
59 
60 template <short_t dim, class T >
61 void gsMaterialMatrixBaseDim<dim,T>::thickness_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
62 {
63  gsMapData<T> map;
64  map.flags = NEED_VALUE;
65  map.points = u;
66  static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
67  m_thickness->piece(patch).eval_into(map.values[0], result);
68 }
69 
70 // Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) covariant basis
71 template <short_t dim, class T>
72 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_spec2cov(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
73 {
74  gsMatrix<T> result(9, u.cols());
75  gsMatrix<T> covbasis,sbasis;
76  gsMatrix<T> tmp = eval3D_pstretchDir(patch,u,z);
77  index_t colIdx;
78  for (index_t k=0; k!=u.cols(); k++)
79  {
80  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
81  {
82  colIdx = j*u.cols()+k;
83  _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
84  sbasis = tmp.reshapeCol(colIdx,3,3);
85  covbasis = m_data.mine().m_gcov_ori;
86  result.col(colIdx) = this->_transformation(covbasis,sbasis).reshape(9,1);
87  }
88  }
89  return result;
90 }
91 
92 // Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) convariant basis
93 template <short_t dim, class T>
94 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_spec2con(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
95 {
96  gsMatrix<T> result(9, u.cols());
97  gsMatrix<T> conbasis,sbasis;
98  gsMatrix<T> tmp = eval3D_pstretchDir(patch,u,z);
99  index_t colIdx;
100  for (index_t k=0; k!=u.cols(); k++)
101  {
102  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
103  {
104  colIdx = j*u.cols()+k;
105  _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
106  sbasis = tmp.reshapeCol(colIdx,3,3);
107  conbasis = m_data.mine().m_gcov_ori;
108  result.col(colIdx) = this->_transformation(conbasis,sbasis).reshape(9,1);
109  }
110  }
111  return result;
112 }
113 
114 // Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) covariant basis
115 template <short_t dim, class T>
116 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_cov2cart(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
117 {
118  gsMatrix<T> result(9, u.cols());
119  gsMatrix<T> tmp, covbasis,cartbasis(3,3);
120  cartbasis.setIdentity();
121  index_t colIdx;
122  for (index_t k=0; k!=u.cols(); k++)
123  {
124  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
125  {
126  colIdx = j*u.cols()+k;
127  _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
128  covbasis = m_data.mine().m_gcov_ori;
129  result.col(colIdx) = this->_transformation(cartbasis,covbasis).reshape(9,1);
130  }
131  }
132  return result;
133 }
134 
135 // Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) convariant basis
136 template <short_t dim, class T>
137 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_con2cart(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
138 {
139  gsMatrix<T> result(9, u.cols());
140  gsMatrix<T> tmp, conbasis,cartbasis(3,3);
141  cartbasis.setIdentity();
142  index_t colIdx;
143  for (index_t k=0; k!=u.cols(); k++)
144  {
145  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
146  {
147  colIdx = j*u.cols()+k;
148  _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
149  conbasis = m_data.mine().m_gcon_ori;
150  result.col(colIdx) = this->_transformation(cartbasis,conbasis).reshape(9,1);
151  }
152  }
153  return result;
154 }
155 
156 template <short_t dim, class T >
157 void gsMaterialMatrixBaseDim<dim,T>::parameters_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
158 {
159  gsMapData<T> map;
160  map.flags = NEED_VALUE;
161  map.points = u;
162  static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
163 
164  gsMatrix<T> tmp;
165  result.resize(m_pars.size(),map.values[0].cols());
166  result.setZero();
167  for (size_t v=0; v!=m_pars.size(); v++)
168  {
169  m_pars[v]->piece(patch).eval_into(map.values[0], tmp);
170  result.row(v) = tmp;
171  }
172 }
173 
174 // Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) convariant basis
175 template <short_t dim, class T >
176 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_deformation(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
177 {
178  this->_computePoints(patch,u);
179 
180  gsMatrix<T> result(9, u.cols() * z.rows());
181  std::pair<gsVector<T>,gsMatrix<T>> res;
182  index_t colIdx;
183  for (index_t k=0; k!=u.cols(); k++)
184  {
185  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
186  {
187  colIdx = j*u.cols()+k;
188  _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
189 
190  gsAsMatrix<T> C = result.reshapeCol(colIdx,3,3);
191  C.setZero();
192  C.block(0,0,2,2) = m_data.mine().m_Gcov_def.block(0,0,2,2);
193  C(2,2) = 1./m_data.mine().m_J0_sq;
194  }
195  }
196  return result;
197 }
198 
199 template <short_t dim, class T>
200 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_pstretch(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
201 {
202  gsMatrix<T> result(3, u.cols() * z.rows());
203  gsMatrix<T> C;
204  std::pair<gsMatrix<T>,gsMatrix<T>> res;
205  gsMatrix<T> tmp = eval3D_deformation(patch,u,z);
206  index_t colIdx;
207  for (index_t k=0; k!= u.cols(); k++)
208  {
209  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
210  {
211  colIdx = j*u.cols()+k;
212  this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
213  C = tmp.reshapeCol(colIdx,3,3);
214  res = this->_evalStretch(C,m_data.mine().m_gcon_ori);
215  result.col(colIdx) = res.first;
216  }
217  }
218  return result;
219 }
220 
221 template <short_t dim, class T>
222 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_pstretchDir(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
223 {
224  gsMatrix<T> result(9, u.cols() * z.rows());
225  gsMatrix<T> C;
226  std::pair<gsMatrix<T>,gsMatrix<T>> res;
227  gsMatrix<T> tmp = eval3D_deformation(patch,u,z);
228  index_t colIdx;
229  for (index_t k=0; k!= u.cols(); k++)
230  {
231  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
232  {
233  colIdx = j*u.cols()+k;
234  this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
235  C = tmp.reshapeCol(colIdx,3,3);
236  res = this->_evalStretch(C,m_data.mine().m_gcon_ori);
237  result.col(colIdx) = res.second.reshape(9,1);
238  }
239  }
240  return result;
241 }
242 
243 template <short_t dim, class T>
244 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_pstrain(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
245 {
246  gsMatrix<T> Emat = eval3D_strain(patch,u,z);
247 
248  gsMatrix<T> result(3, u.cols() * z.rows());
249  result.setZero();
250  gsMatrix<T,3,3> E;
251  std::pair<gsVector<T>,gsMatrix<T>> res;
252  index_t colIdx;
253  for (index_t k=0; k!=u.cols(); k++)
254  {
255  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
256  {
257  colIdx = j * u.cols() + k;
258  E.setZero();
259  E(0,0) = Emat(0,colIdx);
260  E(1,1) = Emat(1,colIdx);
261  E(0,1) = E(1,0) = 0.5*Emat(2,colIdx);
262  E(2,2) = 0;
263  res = this->_evalPStrain(E);
264  result.col(j * u.cols() + k) = res.first;
265  }
266  }
267  return result;
268 }
269 
270 template <short_t dim, class T>
271 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_strain(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
272 {
273  this->_computePoints(patch,u);
274 
275  gsMatrix<T> result(3, u.cols() * z.rows());
276  gsMatrix<T,2,2> E;
277  result.setZero();
278  index_t colIdx;
279  for (index_t k=0; k!=u.cols(); k++)
280  {
281  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
282  {
283  colIdx = j*u.cols()+k;
284  this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
285 
286  E = 0.5 * (m_data.mine().m_Gcov_def.block(0,0,2,2) - m_data.mine().m_Gcov_ori.block(0,0,2,2));
287  result(0,colIdx) = E(0,0);
288  result(1,colIdx) = E(1,1);
289  result(2,colIdx) = E(0,1) + E(1,0);
290  }
291  }
292  return result;
293 }
294 
295 template <short_t dim, class T>
296 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_tensionfield(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
297 {
298  gsMatrix<T> Spmat = this->eval3D_pstress(patch,u,z,MaterialOutput::Generic);
299  gsMatrix<T> Epmat = this->eval3D_pstrain(patch,u,z);
300  gsVector<T> Sp, Ep;
301  gsMatrix<T> result(1, u.cols() * z.rows());
302  index_t colIdx;
303  if (m_options.getInt("TensionField")==0)
304  {
305  for (index_t k=0; k!=u.cols(); k++)
306  {
307  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
308  {
309  colIdx = j*u.cols() + k;
310  Sp = Spmat.col(colIdx);
311  Ep = Epmat.col(colIdx);
312  result.col(colIdx)<<_tensionField<0>(Sp,Ep);
313  }
314  }
315  }
316  else if (m_options.getInt("TensionField")==1)
317  {
318  for (index_t k=0; k!=u.cols(); k++)
319  {
320  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
321  {
322  colIdx = j*u.cols() + k;
323  Sp = Spmat.col(colIdx);
324  Ep = Epmat.col(colIdx);
325  result.col(colIdx)<<_tensionField<1>(Sp,Ep);
326  }
327  }
328  }
329  else if (m_options.getInt("TensionField")==2)
330  {
331  for (index_t k=0; k!=u.cols(); k++)
332  {
333  for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
334  {
335  colIdx = j*u.cols() + k;
336  Sp = Spmat.col(colIdx);
337  Ep = Epmat.col(colIdx);
338  result.col(colIdx)<<_tensionField<2>(Sp,Ep);
339  }
340  }
341  }
342  return result;
343 }
344 
345 template <short_t dim, class T>
346 template <short_t _type>
347 typename std::enable_if<_type==0, index_t>::type
348 gsMaterialMatrixBaseDim<dim,T>::_tensionField(const gsVector<T> & , const gsVector<T> & Ep) const
349 {
350  if (Ep[0] > 0 || (math::abs(Ep[0]) <= Ep.maxCoeff()*1e-12) ) // taut
351  return 1;
352  else if (Ep[1] < 0) // slack
353  return -1;
354  else // wrinkled
355  return 0;
356 }
357 
358 template <short_t dim, class T>
359 template <short_t _type>
360 typename std::enable_if<_type==1, index_t>::type
361 gsMaterialMatrixBaseDim<dim,T>::_tensionField(const gsVector<T> & Sp, const gsVector<T> & ) const
362 {
363  if (Sp[0] > 0 || (math::abs(Sp[0]) <= Sp.maxCoeff()*1e-12) ) // taut
364  return 1;
365  else if (Sp[1] < 0) // slack
366  return -1;
367  else // wrinkled
368  return 0;
369 }
370 
371 template <short_t dim, class T>
372 template <short_t _type>
373 typename std::enable_if<_type==2, index_t>::type
374 gsMaterialMatrixBaseDim<dim,T>::_tensionField(const gsVector<T> & Sp, const gsVector<T> & Ep) const
375 {
376  // See Nakashino 2020
377  // Smin = Sp[0], Smax = Sp[1], S33 = Sp[2]
378  // Emin = Ep[0], Emax = Ep[1], E33 = Ep[2]
379 
380  if (Sp[0] > 0 || (math::abs(Sp[0]) <= Sp.maxCoeff()*1e-12) ) // taut
381  return 1;
382  else if (Ep[1] < 0) // slack
383  return -1;
384  else // wrinkled
385  return 0;
386 }
387 
388 template <short_t dim, class T >
389 void gsMaterialMatrixBaseDim<dim,T>::_computePoints(const index_t patch, const gsMatrix<T> & u) const
390 {
391  gsMatrix<T> tmp;
392 
393  this->_computeMetricUndeformed(patch,u);
394  if (Base::m_defpatches!=nullptr)
395  this->_computeMetricDeformed(patch,u);
396 
397  gsMapData<T> map;
398  map.flags = NEED_VALUE;
399  map.points = u;
400  static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
401 
402  m_thickness->piece(patch).eval_into(map.values[0], m_data.mine().m_Tmat);
403 
404  m_data.mine().m_parmat.resize(m_pars.size(),map.values[0].cols());
405  m_data.mine().m_parmat.setZero();
406 
407  for (size_t v=0; v!=m_pars.size(); v++)
408  {
409  m_pars[v]->piece(patch).eval_into(map.values[0], tmp);
410  m_data.mine().m_parmat.row(v) = tmp;
411  }
412 
413  m_data.mine().m_parvals.resize(m_pars.size());
414 }
415 
416 template <short_t dim, class T>
417 void gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed(const index_t patch, const gsMatrix<T> & u) const
418 {
419  _computeMetricDeformed_impl<dim>(patch,u);
420 }
421 
422 template <short_t dim, class T>
423 template <short_t _dim>
424 typename std::enable_if<_dim==3, void>::type
425 gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed_impl(const index_t patch, const gsMatrix<T> & u) const
426 {
427  gsMapData<T> map;
428  map.flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
429  map.points = u;
430  dynamic_cast<const gsFunction<T>&>(Base::m_defpatches->piece(patch) ).computeMap(map);
431 
432  gsMatrix<T> deriv2(3,3), mixedB(2,2), acov(3,2), acon(3,2), ncov(3,2), Acov(2,2), Acon(2,2), Bcov(2,2);
433  gsMatrix<T> normals;
434  gsVector<T> normal;
435 
436  normals = map.normals;
437  normals.colwise().normalize();
438  m_data.mine().m_normal_def_mat = normals;
439 
440  m_data.mine().m_Acov_def_mat.resize(4,map.points.cols()); m_data.mine().m_Acov_def_mat.setZero();
441  m_data.mine().m_Acon_def_mat.resize(4,map.points.cols()); m_data.mine().m_Acon_def_mat.setZero();
442  m_data.mine().m_Bcov_def_mat.resize(4,map.points.cols()); m_data.mine().m_Bcov_def_mat.setZero();
443 
444  m_data.mine().m_acov_def_mat.resize(2*3,map.points.cols()); m_data.mine().m_acov_def_mat.setZero();
445  m_data.mine().m_acon_def_mat.resize(2*3,map.points.cols()); m_data.mine().m_acon_def_mat.setZero();
446  m_data.mine().m_ncov_def_mat.resize(2*3,map.points.cols()); m_data.mine().m_ncov_def_mat.setZero();
447 
448  for (index_t k=0; k!= map.points.cols(); k++)
449  {
450  acov = map.jacobian(k);
451 
452  Acov = acov.transpose() * acov;
453  Acon = Acov.inverse();
454 
455  // Construct metric tensor b = [d11c*n, d12c*n ; d21c*n, d22c*n]
456  deriv2 = map.deriv2(k).reshaped(3,3);
457  normal = normals.col(k);
458 
459  Bcov(0,0) = deriv2.row(0).dot(normal);
460  Bcov(1,1) = deriv2.row(1).dot(normal);
461  Bcov(0,1) = Bcov(1,0) = deriv2.row(2).dot(normal);
462 
463  // Construct basis
464  for (index_t i=0; i < 2; i++)
465  acon.col(i) = Acon(i,0)*acov.col(0) + Acon(i,1)*acov.col(1);
466 
467  // Mixed tensor
468  for (index_t i=0; i < 2; i++)
469  for (index_t j=0; j < 2; j++)
470  mixedB(i,j) = Acon(i,0)*Bcov(0,j) + Acon(i,1)*Bcov(1,j);
471 
472  for (index_t i=0; i < 2; i++)
473  ncov.col(i) = -mixedB(0,i)*acov.col(0) -mixedB(1,i)*acov.col(1);
474 
475  // Assign members
476  m_data.mine().m_acov_def_mat.reshapeCol(k,3,2) = acov;
477  m_data.mine().m_acon_def_mat.reshapeCol(k,3,2) = acon;
478  m_data.mine().m_ncov_def_mat.reshapeCol(k,3,2) = ncov;
479  m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2) = Acov;
480  m_data.mine().m_Acon_def_mat.reshapeCol(k,2,2) = Acon;
481  m_data.mine().m_Bcov_def_mat.reshapeCol(k,2,2) = Bcov;
482  }
483 }
484 
485 template <short_t dim, class T>
486 template <short_t _dim>
487 typename std::enable_if<_dim==2, void>::type
488 gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed_impl(const index_t patch, const gsMatrix<T> & u) const
489 {
490  gsMapData<T> map;
491  map.flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
492  map.points = u;
493  dynamic_cast<const gsFunction<T>&>(Base::m_defpatches->piece(patch) ).computeMap(map);
494 
495  gsMatrix<T> tmp;
496  gsMatrix<T> acov(2,2), acon(2,2), Acov(2,2), Acon(2,2);
497 
499  m_data.mine().m_Acov_def_mat.resize(4,map.points.cols()); m_data.mine().m_Acov_def_mat.setZero();
500  m_data.mine().m_Acon_def_mat.resize(4,map.points.cols()); m_data.mine().m_Acon_def_mat.setZero();
501  m_data.mine().m_Bcov_def_mat.resize(4,map.points.cols()); m_data.mine().m_Bcov_def_mat.setZero();
502 
503  m_data.mine().m_acov_def_mat.resize(2*2,map.points.cols()); m_data.mine().m_acov_def_mat.setZero();
504  m_data.mine().m_acon_def_mat.resize(2*2,map.points.cols()); m_data.mine().m_acon_def_mat.setZero();
505  m_data.mine().m_ncov_def_mat.resize(2*2,map.points.cols()); m_data.mine().m_ncov_def_mat.setZero();
506 
507  gsMatrix<T> zero(2,2); zero.setZero();
508  for (index_t k=0; k!= map.points.cols(); k++)
509  {
510  acov = map.jacobian(k);
511 
512  Acov = acov.transpose() * acov;
513  Acon = Acov.inverse();
514 
515  // Construct basis
516  for (index_t i=0; i < 2; i++)
517  acon.col(i) = Acon(i,0)*acov.col(0) + Acon(i,1)*acov.col(1);
518 
519  // Assign members
520  m_data.mine().m_acov_def_mat.reshapeCol(k,2,2) = acov;
521  m_data.mine().m_acon_def_mat.reshapeCol(k,2,2) = acon;
522  m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2) = Acov;
523  m_data.mine().m_Acon_def_mat.reshapeCol(k,2,2) = Acon;
524 
525  // Since setZero above does not work
526  m_data.mine().m_Bcov_def_mat.reshapeCol(k,2,2) = zero;
527  m_data.mine().m_ncov_def_mat.reshapeCol(k,2,2) = zero;
528  }
529 }
530 
531 
532 //--------------------------------------------------------------------------------------------------------------------------------------
533 
534 template <short_t dim, class T>
535 void gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed(const index_t patch, const gsMatrix<T> & u) const
536 {
537  _computeMetricUndeformed_impl<dim>(patch,u);
538 }
539 
540 template <short_t dim, class T>
541 template <short_t _dim>
542 typename std::enable_if<_dim==3, void>::type
543 gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed_impl(const index_t patch, const gsMatrix<T> & u) const
544 {
545  gsMapData<T> map;
546  map.flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
547  map.points = u;
548  dynamic_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
549 
550  gsMatrix<T> deriv2(3,3), mixedB(2,2), acov(3,2), acon(3,2), ncov(3,2), Acov(2,2), Acon(2,2), Bcov(2,2);
551  gsMatrix<T> normals;
552  gsVector<T> normal;
553 
554  normals = map.normals;
555  normals.colwise().normalize();
556  m_data.mine().m_normal_ori_mat = normals;
557 
558  m_data.mine().m_Acov_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Acov_ori_mat.setZero();
559  m_data.mine().m_Acon_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Acon_ori_mat.setZero();
560  m_data.mine().m_Bcov_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Bcov_ori_mat.setZero();
561 
562  m_data.mine().m_acov_ori_mat.resize(2*3,map.points.cols()); m_data.mine().m_acov_ori_mat.setZero();
563  m_data.mine().m_acon_ori_mat.resize(2*3,map.points.cols()); m_data.mine().m_acon_ori_mat.setZero();
564  m_data.mine().m_ncov_ori_mat.resize(2*3,map.points.cols()); m_data.mine().m_ncov_ori_mat.setZero();
565 
566  for (index_t k=0; k!= map.points.cols(); k++)
567  {
568  acov = map.jacobian(k);
569 
570  Acov = acov.transpose() * acov;
571  Acon = Acov.inverse();
572 
573  // Construct metric tensor b = [d11c*n, d12c*n ; d21c*n, d22c*n]
574  deriv2 = map.deriv2(k).reshaped(3,3);
575  normal = normals.col(k);
576 
577  Bcov(0,0) = deriv2.row(0).dot(normal);
578  Bcov(1,1) = deriv2.row(1).dot(normal);
579  Bcov(0,1) = Bcov(1,0) = deriv2.row(2).dot(normal);
580 
581  // Construct basis
582  for (index_t i=0; i < 2; i++)
583  acon.col(i) = Acon(i,0)*acov.col(0) + Acon(i,1)*acov.col(1);
584 
585  // Mixed tensor
586  for (index_t i=0; i < 2; i++)
587  for (index_t j=0; j < 2; j++)
588  mixedB(i,j) = Acon(i,0)*Bcov(0,j) + Acon(i,1)*Bcov(1,j);
589 
590  for (index_t i=0; i < 2; i++)
591  ncov.col(i) = -mixedB(0,i)*acov.col(0) -mixedB(1,i)*acov.col(1);
592 
593  // Assign members
594  m_data.mine().m_acov_ori_mat.reshapeCol(k,3,2) = acov;
595  m_data.mine().m_acon_ori_mat.reshapeCol(k,3,2) = acon;
596  m_data.mine().m_ncov_ori_mat.reshapeCol(k,3,2) = ncov;
597  m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2) = Acov;
598  m_data.mine().m_Acon_ori_mat.reshapeCol(k,2,2) = Acon;
599  m_data.mine().m_Bcov_ori_mat.reshapeCol(k,2,2) = Bcov;
600  }
601 }
602 
603 template <short_t dim, class T>
604 template <short_t _dim>
605 typename std::enable_if<_dim==2, void>::type
606 gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed_impl(const index_t patch, const gsMatrix<T> & u) const
607 {
608  gsMapData<T> map;
609  map.flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
610  map.points = u;
611  dynamic_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
612 
613  gsMatrix<T> tmp;
614  gsMatrix<T> acov(2,2), acon(2,2), Acov(2,2), Acon(2,2);
615 
616  m_data.mine().m_Acov_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Acov_ori_mat.setZero();
617  m_data.mine().m_Acon_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Acon_ori_mat.setZero();
618  m_data.mine().m_Bcov_ori_mat.resize(4,map.points.cols()); m_data.mine().m_Bcov_ori_mat.setZero();
619 
620  m_data.mine().m_acov_ori_mat.resize(2*2,map.points.cols()); m_data.mine().m_acov_ori_mat.setZero();
621  m_data.mine().m_acon_ori_mat.resize(2*2,map.points.cols()); m_data.mine().m_acon_ori_mat.setZero();
622  m_data.mine().m_ncov_ori_mat.resize(2*2,map.points.cols()); m_data.mine().m_ncov_ori_mat.setZero();
623 
624  gsMatrix<T> zero(2,2); zero.setZero();
625  for (index_t k=0; k!= map.points.cols(); k++)
626  {
627  acov = map.jacobian(k);
628 
629  Acov = acov.transpose() * acov;
630  Acon = Acov.inverse();
631 
632  // Construct basis
633  for (index_t i=0; i < 2; i++)
634  acon.col(i) = Acon(i,0)*acov.col(0) + Acon(i,1)*acov.col(1);
635 
636  // Assign members
637  m_data.mine().m_acov_ori_mat.reshapeCol(k,2,2) = acov;
638  m_data.mine().m_acon_ori_mat.reshapeCol(k,2,2) = acon;
639  m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2) = Acov;
640  m_data.mine().m_Acon_ori_mat.reshapeCol(k,2,2) = Acon;
641 
642  // Since setZero above does not work
643  m_data.mine().m_Bcov_ori_mat.reshapeCol(k,2,2) = zero;
644  m_data.mine().m_ncov_ori_mat.reshapeCol(k,2,2) = zero;
645  }
646 }
647 
648 //--------------------------------------------------------------------------------------------------------------------------------------
649 
650 template <short_t dim, class T>
651 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getAcov_def(index_t k, T z) const
652 {
653  GISMO_ENSURE(m_data.mine().m_Acov_def_mat.cols()!=0,"Is the metric initialized?");
654  return m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2);
655 }
656 
657 template <short_t dim, class T>
658 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getAcon_def(index_t k, T z) const
659 {
660  GISMO_ENSURE(m_data.mine().m_Acon_def_mat.cols()!=0,"Is the metric initialized?");
661  return m_data.mine().m_Acon_def_mat.reshapeCol(k,2,2);
662 }
663 
664 template <short_t dim, class T>
665 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getBcov_def(index_t k, T z) const
666 {
667  return _getBcov_def_impl<dim>(k,z);
668 }
669 
670 template <short_t dim, class T>
671 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getncov_def(index_t k, T z) const
672 {
673  return _getncov_def_impl<dim>(k,z);
674 }
675 
676 template <short_t dim, class T>
677 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcov_def(index_t k, T z) const
678 {
679  return _getGcov_def_impl<dim>(k,z);
680 }
681 
682 template <short_t dim, class T>
683 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcon_def(index_t k, T z) const
684 {
685  gsMatrix<T> Gcov_def = _getGcov_def(k,z);
686  return Gcov_def.inverse();
687 }
688 
689 template <short_t dim, class T>
690 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getacov_def(index_t k, T z) const
691 {
692  GISMO_ENSURE(m_data.mine().m_acov_def_mat.cols()!=0,"Is the basis initialized?");
693  return m_data.mine().m_acov_def_mat.reshapeCol(k,3,2);
694 }
695 
696 template <short_t dim, class T>
697 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getacon_def(index_t k, T z) const
698 {
699  GISMO_ENSURE(m_data.mine().m_acon_def_mat.cols()!=0,"Is the basis initialized?");
700  return m_data.mine().m_acon_def_mat.reshapeCol(k,3,2);
701 }
702 
703 template <short_t dim, class T>
704 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcov_def(index_t k, T z) const
705 {
706  return _getgcov_def_impl<dim>(k,z);
707 }
708 
709 template <short_t dim, class T>
710 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcon_def(index_t k, T z) const
711 {
712  gsMatrix<T> gcov_def = _getgcov_def(k,z);
713  gsMatrix<T> Gcon_def = _getGcon_def(k,z);
714  gsMatrix<T> gcon_def(3,3);
715  for (index_t c = 0; c!=3; c++)
716  {
717  gcon_def.col(c) = Gcon_def(c,0) * gcov_def.col(0)
718  + Gcon_def(c,1) * gcov_def.col(1)
719  + Gcon_def(c,2) * gcov_def.col(2);
720  }
721  return gcon_def;
722 }
723 
724 template <short_t dim, class T>
725 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getAcov_ori(index_t k, T z) const
726 {
727  GISMO_ENSURE(m_data.mine().m_Acov_ori_mat.cols()!=0,"Is the metric initialized?");
728  return m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2);
729 }
730 
731 template <short_t dim, class T>
732 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getAcon_ori(index_t k, T z) const
733 {
734  GISMO_ENSURE(m_data.mine().m_Acon_ori_mat.cols()!=0,"Is the metric initialized?");
735  return m_data.mine().m_Acon_ori_mat.reshapeCol(k,2,2);
736 }
737 
738 template <short_t dim, class T>
739 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getBcov_ori(index_t k, T z) const
740 {
741  return _getBcov_ori_impl<dim>(k,z);
742 }
743 
744 template <short_t dim, class T>
745 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getncov_ori(index_t k, T z) const
746 {
747  return _getncov_ori_impl<dim>(k,z);
748 }
749 
750 template <short_t dim, class T>
751 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcov_ori(index_t k, T z) const
752 {
753  return _getGcov_ori_impl<dim>(k,z);
754 }
755 
756 template <short_t dim, class T>
757 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcon_ori(index_t k, T z) const
758 {
759  gsMatrix<T> Gcov_ori = _getGcov_ori(k,z);
760  return Gcov_ori.inverse();
761 }
762 
763 template <short_t dim, class T>
764 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getacov_ori(index_t k, T z) const
765 {
766  GISMO_ENSURE(m_data.mine().m_acov_ori_mat.cols()!=0,"Is the basis initialized?");
767  return m_data.mine().m_acov_ori_mat.reshapeCol(k,3,2);
768 }
769 
770 template <short_t dim, class T>
771 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getacon_ori(index_t k, T z) const
772 {
773  GISMO_ENSURE(m_data.mine().m_acon_ori_mat.cols()!=0,"Is the basis initialized?");
774  return m_data.mine().m_acon_ori_mat.reshapeCol(k,3,2);
775 }
776 
777 template <short_t dim, class T>
778 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcov_ori(index_t k, T z) const
779 {
780  return _getgcov_ori_impl<dim>(k,z);
781 }
782 
783 template <short_t dim, class T>
784 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcon_ori(index_t k, T z) const
785 {
786  gsMatrix<T> gcov_ori = _getgcov_ori(k,z);
787  gsMatrix<T> Gcon_ori = _getGcon_ori(k,z);
788  gsMatrix<T> gcon_ori(3,3);
789  for (index_t c = 0; c!=3; c++)
790  {
791  gcon_ori.col(c) = Gcon_ori(c,0) * gcov_ori.col(0)
792  + Gcon_ori(c,1) * gcov_ori.col(1)
793  + Gcon_ori(c,2) * gcov_ori.col(2);
794  }
795  return gcon_ori;
796 }
797 
798 //--------------------------------------------------------------------------------------------------------------------------------------
799 
800 template <short_t dim, class T>
801 template <short_t _dim>
802 typename std::enable_if<_dim==2, gsMatrix<T>>::type
803 gsMaterialMatrixBaseDim<dim,T>::_getBcov_def_impl(index_t k, T z) const
804 {
805  GISMO_NO_IMPLEMENTATION;
806 }
807 
808 template <short_t dim, class T>
809 template <short_t _dim>
810 typename std::enable_if<_dim==3, gsMatrix<T>>::type
811 gsMaterialMatrixBaseDim<dim,T>::_getBcov_def_impl(index_t k, T z) const
812 {
813  GISMO_ENSURE(m_data.mine().m_Bcov_def_mat.cols()!=0,"Is the metric initialized?");
814  return m_data.mine().m_Bcov_def_mat.reshapeCol(k,2,2);
815 }
816 
817 template <short_t dim, class T>
818 template <short_t _dim>
819 typename std::enable_if<_dim==2, gsMatrix<T>>::type
820 gsMaterialMatrixBaseDim<dim,T>::_getncov_def_impl(index_t k, T z) const
821 {
822  GISMO_NO_IMPLEMENTATION;
823 }
824 
825 template <short_t dim, class T>
826 template <short_t _dim>
827 typename std::enable_if<_dim==3, gsMatrix<T>>::type
828 gsMaterialMatrixBaseDim<dim,T>::_getncov_def_impl(index_t k, T z) const
829 {
830  GISMO_ENSURE(m_data.mine().m_ncov_def_mat.cols()!=0,"Is the basis initialized?");
831  return m_data.mine().m_ncov_def_mat.reshapeCol(k,3,2);
832 }
833 
834 template <short_t dim, class T>
835 template <short_t _dim>
836 typename std::enable_if<_dim==2, gsMatrix<T>>::type
837 gsMaterialMatrixBaseDim<dim,T>::_getGcov_def_impl(index_t k, T z) const
838 {
839  // metrics
840  gsMatrix<T> Acov_def = _getAcov_def(k,z);
841 
842  // Compute full metric
843  gsMatrix<T> Gcov_def(3,3);
844  Gcov_def.setZero();
845  Gcov_def.block(0,0,2,2)= Acov_def;
846  Gcov_def(2,2) = 1.0;
847  return Gcov_def;
848 }
849 
850 template <short_t dim, class T>
851 template <short_t _dim>
852 typename std::enable_if<_dim==3, gsMatrix<T>>::type
853 gsMaterialMatrixBaseDim<dim,T>::_getGcov_def_impl(index_t k, T z) const
854 {
855  // metrics
856  gsMatrix<T> Acov_def = _getAcov_def(k,z);
857  gsMatrix<T> Bcov_def = _getBcov_def(k,z);
858  gsMatrix<T> ncov_def = _getncov_def(k,z);
859 
860  // Compute full metric
861  gsMatrix<T> Gcov_def(3,3);
862  Gcov_def.setZero();
863  Gcov_def.block(0,0,2,2)= Acov_def - 2.0 * z * Bcov_def + z*z * ncov_def.transpose()*ncov_def;
864  Gcov_def(2,2) = 1.0;
865  return Gcov_def;
866 }
867 
868 template <short_t dim, class T>
869 template <short_t _dim>
870 typename std::enable_if<_dim==2, gsMatrix<T>>::type
871 gsMaterialMatrixBaseDim<dim,T>::_getgcov_def_impl(index_t k, T z) const
872 {
873  gsMatrix<T> normal(3,1);
874  normal << 0,0,1;
875  gsMatrix<T> gcov_def(3,3);
876  gcov_def.setZero();
877  gcov_def.block(0,0,2,2) = _getacov_def(k,z);
878  gcov_def.col(2) = normal;
879  return gcov_def;
880 }
881 
882 template <short_t dim, class T>
883 template <short_t _dim>
884 typename std::enable_if<_dim==3, gsMatrix<T>>::type
885 gsMaterialMatrixBaseDim<dim,T>::_getgcov_def_impl(index_t k, T z) const
886 {
887  GISMO_ENSURE(m_data.mine().m_normal_def_mat.cols()!=0,"Is the basis initialized?");
888  gsMatrix<T> normal = m_data.mine().m_normal_def_mat.reshapeCol(k,3,1);
889  gsMatrix<T> acov_def = _getacov_def(k,z);
890  gsMatrix<T> ncov_def = _getncov_def(k,z);
891  gsMatrix<T> gcov_def(3,3);
892  gcov_def.setZero();
893  gcov_def.leftCols(2) = acov_def + z * ncov_def;
894  gcov_def.col(2) = normal;
895  return gcov_def;
896 }
897 
898 template <short_t dim, class T>
899 template <short_t _dim>
900 typename std::enable_if<_dim==2, gsMatrix<T>>::type
901 gsMaterialMatrixBaseDim<dim,T>::_getBcov_ori_impl(index_t k, T z) const
902 {
903  GISMO_NO_IMPLEMENTATION;
904 }
905 
906 template <short_t dim, class T>
907 template <short_t _dim>
908 typename std::enable_if<_dim==3, gsMatrix<T>>::type
909 gsMaterialMatrixBaseDim<dim,T>::_getBcov_ori_impl(index_t k, T z) const
910 {
911  GISMO_ENSURE(m_data.mine().m_Bcov_ori_mat.cols()!=0,"Is the metric initialized?");
912  return m_data.mine().m_Bcov_ori_mat.reshapeCol(k,2,2);
913 }
914 
915 template <short_t dim, class T>
916 template <short_t _dim>
917 typename std::enable_if<_dim==2, gsMatrix<T>>::type
918 gsMaterialMatrixBaseDim<dim,T>::_getncov_ori_impl(index_t k, T z) const
919 {
920  GISMO_NO_IMPLEMENTATION;
921 }
922 
923 template <short_t dim, class T>
924 template <short_t _dim>
925 typename std::enable_if<_dim==3, gsMatrix<T>>::type
926 gsMaterialMatrixBaseDim<dim,T>::_getncov_ori_impl(index_t k, T z) const
927 {
928  GISMO_ENSURE(m_data.mine().m_ncov_ori_mat.cols()!=0,"Is the basis initialized?");
929  return m_data.mine().m_ncov_ori_mat.reshapeCol(k,3,2);
930 }
931 
932 template <short_t dim, class T>
933 template <short_t _dim>
934 typename std::enable_if<_dim==2, gsMatrix<T>>::type
935 gsMaterialMatrixBaseDim<dim,T>::_getGcov_ori_impl(index_t k, T z) const
936 {
937  // metrics
938  gsMatrix<T> Acov_ori = _getAcov_ori(k,z);
939 
940  // Compute full metric
941  gsMatrix<T> Gcov_ori(3,3);
942  Gcov_ori.setZero();
943  Gcov_ori.block(0,0,2,2)= Acov_ori;
944  Gcov_ori(2,2) = 1.0;
945  return Gcov_ori;
946 }
947 
948 template <short_t dim, class T>
949 template <short_t _dim>
950 typename std::enable_if<_dim==3, gsMatrix<T>>::type
951 gsMaterialMatrixBaseDim<dim,T>::_getGcov_ori_impl(index_t k, T z) const
952 {
953  // metrics
954  gsMatrix<T> Acov_ori = _getAcov_ori(k,z);
955  gsMatrix<T> Bcov_ori = _getBcov_ori(k,z);
956  gsMatrix<T> ncov_ori = _getncov_ori(k,z);
957 
958  // Compute full metric
959  gsMatrix<T> Gcov_ori(3,3);
960  Gcov_ori.setZero();
961  Gcov_ori.block(0,0,2,2)= Acov_ori - 2.0 * z * Bcov_ori + z*z * ncov_ori.transpose()*ncov_ori;
962  Gcov_ori(2,2) = 1.0;
963  return Gcov_ori;
964 }
965 
966 template <short_t dim, class T>
967 template <short_t _dim>
968 typename std::enable_if<_dim==2, gsMatrix<T>>::type
969 gsMaterialMatrixBaseDim<dim,T>::_getgcov_ori_impl(index_t k, T z) const
970 {
971  gsMatrix<T> normal(3,1);
972  normal << 0,0,1;
973  gsMatrix<T> gcov_ori(3,3);
974  gcov_ori.setZero();
975  gcov_ori.block(0,0,2,2) = _getacov_ori(k,z);
976  gcov_ori.col(2) = normal;
977  return gcov_ori;
978 }
979 
980 template <short_t dim, class T>
981 template <short_t _dim>
982 typename std::enable_if<_dim==3, gsMatrix<T>>::type
983 gsMaterialMatrixBaseDim<dim,T>::_getgcov_ori_impl(index_t k, T z) const
984 {
985  GISMO_ENSURE(m_data.mine().m_normal_ori_mat.cols()!=0,"Is the basis initialized?");
986  gsMatrix<T> normal = m_data.mine().m_normal_ori_mat.reshapeCol(k,3,1);
987  gsMatrix<T> acov_ori = _getacov_ori(k,z);
988  gsMatrix<T> ncov_ori = _getncov_ori(k,z);
989  gsMatrix<T> gcov_ori(3,3);
990  gcov_ori.setZero();
991  gcov_ori.leftCols(2) = acov_ori + z * ncov_ori;
992  gcov_ori.col(2) = normal;
993  return gcov_ori;
994 }
995 
996 //--------------------------------------------------------------------------------------------------------------------------------------
997 
998 template <short_t dim, class T>
999 void gsMaterialMatrixBaseDim<dim,T>::_getMetric(index_t k, T z, const gsMatrix<T> & C) const
1000 {
1001  this->_getMetricDeformed(C);
1002  this->_getMetricUndeformed(k,z);
1003 
1004  T ratio;
1005  T det_ori = m_data.mine().m_Gcov_ori.determinant();
1006  T det_def = m_data.mine().m_Gcov_def.determinant();
1007 
1008  if ((det_ori==0 && det_def==0) || (math::isnan(det_ori) && math::isnan(det_def)))
1009  {
1010  gsWarn<<"Jacobian determinant is undefined: J^2 = det(Gcov_def) / det(Gcov_ori) = "<<det_def<<"/"<<det_ori<<"! J^2 is set to 1";
1011  ratio = 1;
1012  }
1013  else
1014  ratio = det_def / det_ori;
1015 
1016  GISMO_ENSURE(ratio >= 0, "Jacobian determinant is negative! det(Gcov_def) = "<<det_def<<"; det(Gcov_ori) = "<<det_ori);
1017  m_data.mine().m_J0_sq = ratio;
1018 }
1019 
1020 template <short_t dim, class T>
1021 void gsMaterialMatrixBaseDim<dim,T>::_getMetric(const index_t k, const T z) const
1022 {
1023  this->_getMetricDeformed(k,z);
1024  this->_getMetricUndeformed(k,z);
1025 
1026  T ratio;
1027  T det_ori = m_data.mine().m_Gcov_ori.determinant();
1028  T det_def = m_data.mine().m_Gcov_def.determinant();
1029 
1030  if ((det_ori==0 && det_def==0) || (math::isnan(det_ori) && math::isnan(det_def)))
1031  {
1032  gsWarn<<"Jacobian determinant is undefined: J^2 = det(Gcov_def) / det(Gcov_ori) = "<<det_def<<"/"<<det_ori<<"! J^2 is set to 1";
1033  ratio = 1;
1034  }
1035  else
1036  ratio = det_def / det_ori;
1037 
1038  GISMO_ENSURE(ratio >= 0, "Jacobian determinant is negative! det(Gcov_def) = "<<det_def<<"; det(Gcov_ori) = "<<det_ori<<"\nGcov_def = "<<m_data.mine().m_Gcov_def<<"\n"<<"Acov_def = "<<m_data.mine().m_Acov_def<<"\nBcov_def = "<<m_data.mine().m_Bcov_def);
1039  m_data.mine().m_J0_sq = ratio;
1040 }
1041 
1042 //--------------------------------------------------------------------------------------------------------------------------------------
1043 
1044 template <short_t dim, class T>
1045 void gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed(const gsMatrix<T> & C) const
1046 {
1047  // Compute full metric
1048  gsMatrix<T> Gcov_def(3,3), Gcon_def(3,3);
1049  Gcov_def.setZero();
1050  Gcon_def.setZero();
1051  Gcov_def(0,0) = C(0,0);
1052  Gcov_def(1,1) = C(1,0);
1053  Gcov_def(0,1) =
1054  Gcov_def(1,0) = C(2,0);
1055  Gcov_def(2,2) =
1056  Gcon_def(2,2) = 1.0;
1057  Gcon_def.block(0,0,2,2) = Gcov_def.block(0,0,2,2).inverse();
1058 
1059  m_data.mine().m_Gcov_def = Gcov_def;
1060  m_data.mine().m_Gcon_def = Gcon_def;
1061  // CLear other members on the deformed geometry
1062  m_data.mine().m_Acov_def.setZero();
1063  m_data.mine().m_Acon_def.setZero();
1064  m_data.mine().m_Bcov_def.setZero();
1065  m_data.mine().m_Bcon_def.setZero();
1066  m_data.mine().m_acov_def.setZero();
1067  m_data.mine().m_acon_def.setZero();
1068  m_data.mine().m_gcov_def.setZero();
1069  m_data.mine().m_gcon_def.setZero();
1070 }
1071 
1072 template <short_t dim, class T>
1073 void gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed(const index_t k, const T z) const
1074 {
1075  _getMetricDeformed_impl<dim>(k,z);
1076 }
1077 
1078 template <short_t dim, class T>
1079 template <short_t _dim>
1080 typename std::enable_if<_dim==3, void>::type
1081 gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed_impl(const index_t k, const T z) const
1082 {
1083  GISMO_ENSURE(m_data.mine().m_Acov_def_mat.cols()!=0,"Is the metric initialized?");
1084  GISMO_ENSURE(m_data.mine().m_Acon_def_mat.cols()!=0,"Is the metric initialized?");
1085  GISMO_ENSURE(m_data.mine().m_Bcov_def_mat.cols()!=0,"Is the metric initialized?");
1086  GISMO_ENSURE(m_data.mine().m_ncov_def_mat.cols()!=0,"Is the basis initialized?");
1087  GISMO_ENSURE(m_data.mine().m_acov_def_mat.cols()!=0,"Is the basis initialized?");
1088  GISMO_ENSURE(m_data.mine().m_acon_def_mat.cols()!=0,"Is the basis initialized?");
1089  GISMO_ENSURE(m_data.mine().m_normal_def_mat.cols()!=0,"Is the basis initialized?");
1090 
1091  gsMatrix<T> Acov_def, Acon_def, Bcov_def, ncov_def, Gcov_def(3,3), Gcon_def(3,3),
1092  acov_def, acon_def, normal(3,1), gcov_def(3,3), gcon_def(3,3);
1093  // Get metric information
1094  Acov_def = m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2);
1095  Acon_def = m_data.mine().m_Acon_def_mat.reshapeCol(k,2,2);
1096  Bcov_def = m_data.mine().m_Bcov_def_mat.reshapeCol(k,2,2);
1097  ncov_def = m_data.mine().m_ncov_def_mat.reshapeCol(k,3,2);
1098 
1099  // Compute full metric
1100  Gcov_def.setZero();
1101  Gcov_def.block(0,0,2,2)= Acov_def - 2.0 * z * Bcov_def + z*z * ncov_def.transpose()*ncov_def;
1102  Gcov_def(2,2) = 1.0;
1103  Gcon_def = Gcov_def.inverse();
1104 
1105  // Assign members
1106  m_data.mine().m_Acov_def = Acov_def;
1107  m_data.mine().m_Acon_def = Acon_def;
1108  m_data.mine().m_Bcov_def = Bcov_def;
1109  m_data.mine().m_ncov_def = ncov_def;
1110  m_data.mine().m_Gcov_def = Gcov_def;
1111  m_data.mine().m_Gcon_def = Gcon_def;
1112 
1113  // Get basis vectors
1114  acov_def = m_data.mine().m_acov_def_mat.reshapeCol(k,3,2);
1115  acon_def = m_data.mine().m_acon_def_mat.reshapeCol(k,3,2);
1116  normal = m_data.mine().m_normal_def_mat.reshapeCol(k,3,1);
1117 
1118  // Compute g_cov
1119  gcov_def.setZero();
1120  gcov_def.leftCols(2) = acov_def + z * ncov_def;
1121  gcov_def.col(2) = normal;
1122 
1123  // Compute g_con
1124  for (index_t c = 0; c!=3; c++)
1125  gcon_def.col(c) = Gcon_def(c,0) * gcov_def.col(0) + Gcon_def(c,1) * gcov_def.col(1) + Gcon_def(c,2) * gcov_def.col(2);
1126 
1127  // Assign members
1128  m_data.mine().m_acov_def = acov_def;
1129  m_data.mine().m_acon_def = acov_def;
1130  m_data.mine().m_gcov_def = gcov_def;
1131  m_data.mine().m_gcon_def = gcon_def;
1132 
1133  // // create a Linearised covariant tensor for SvK models
1134  // m_data.mine().m_Gcov_def_L = m_data.mine().m_Gcov_def;
1135  // m_data.mine().m_Gcov_def.block(0,0,2,2) -= z*z * m_data.mine().m_ncov_def.transpose()*m_data.mine().m_ncov_def;
1136 }
1137 
1138 template <short_t dim, class T>
1139 template <short_t _dim>
1140 typename std::enable_if<_dim==2, void>::type
1141 gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed_impl(const index_t k, const T z) const
1142 {
1143  GISMO_ENSURE(m_data.mine().m_Acov_def_mat.cols()!=0,"Is the metric initialized?");
1144  GISMO_ENSURE(m_data.mine().m_Acon_def_mat.cols()!=0,"Is the metric initialized?");
1145  GISMO_ENSURE(m_data.mine().m_Bcov_def_mat.cols()!=0,"Is the metric initialized?");
1146  GISMO_ENSURE(m_data.mine().m_ncov_def_mat.cols()!=0,"Is the basis initialized?");
1147  GISMO_ENSURE(m_data.mine().m_acov_def_mat.cols()!=0,"Is the basis initialized?");
1148  GISMO_ENSURE(m_data.mine().m_acon_def_mat.cols()!=0,"Is the basis initialized?");
1149 
1150  GISMO_ENSURE(m_data.mine().m_acov_def_mat.cols()!=0,"Is the basis initialized?");
1151  GISMO_ENSURE(m_data.mine().m_acon_def_mat.cols()!=0,"Is the basis initialized?");
1152 
1153  gsMatrix<T> Acov_def, Acon_def, Gcov_def(3,3), Gcon_def(3,3),
1154  acov_def, acon_def, normal(3,1), gcov_def(3,3), gcon_def(3,3);
1155 
1156  // Get metric information
1157  Acov_def = m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2);
1158  Acon_def = m_data.mine().m_Acov_def_mat.reshapeCol(k,2,2);
1159 
1160  // Compute full metric
1161  Gcov_def.setZero();
1162  Gcov_def.block(0,0,2,2)= Acov_def;;
1163  Gcov_def(2,2) = 1.0;
1164  Gcon_def = Gcov_def.inverse();
1165 
1166  // Assign members
1167  m_data.mine().m_Acov_def = Acov_def;
1168  m_data.mine().m_Acon_def = Acon_def;
1169  m_data.mine().m_Bcov_def = m_data.mine().m_Bcov_def_mat.reshapeCol(k,2,2);
1170  m_data.mine().m_ncov_def = m_data.mine().m_ncov_def_mat.reshapeCol(k,2,2);
1171  m_data.mine().m_Gcov_def = Gcov_def;
1172  m_data.mine().m_Gcon_def = Gcon_def;
1173 
1174  // Get basis vectors
1175  acov_def = m_data.mine().m_acov_def_mat.reshapeCol(k,2,2);
1176  acon_def = m_data.mine().m_acon_def_mat.reshapeCol(k,2,2);
1177  normal << 0,0,1;
1178 
1179  // Compute g_cov
1180  gcov_def.setZero();
1181  gcov_def.block(0,0,2,2) = acov_def;
1182  gcov_def.col(2) = normal;
1183 
1184  // Compute g_con
1185  for (index_t c = 0; c!=3; c++)
1186  gcon_def.col(c) = Gcon_def(c,0) * gcov_def.col(0) + Gcon_def(c,1) * gcov_def.col(1) + Gcon_def(c,2) * gcov_def.col(2);
1187 
1188  // Assign members
1189  m_data.mine().m_acov_def = acov_def;
1190  m_data.mine().m_acon_def = acon_def;
1191  m_data.mine().m_gcov_def = gcov_def;
1192  m_data.mine().m_gcon_def = gcon_def;
1193 
1194  // // create a Linearised covariant tensor for SvK models
1195  // m_data.mine().m_Gcov_def_L = m_data.mine().m_Gcov_def;
1196  // m_data.mine().m_Gcov_def.block(0,0,2,2) -= z*z * m_data.mine().m_ncov_def.transpose()*m_data.mine().m_ncov_def;
1197 }
1198 
1199 //--------------------------------------------------------------------------------------------------------------------------------------
1200 
1201 template <short_t dim, class T>
1202 void gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed(const index_t k, const T z) const
1203 {
1204  _getMetricUndeformed_impl<dim>(k,z);
1205 }
1206 
1207 template <short_t dim, class T>
1208 template <short_t _dim>
1209 typename std::enable_if<_dim==3, void>::type
1210 gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed_impl(const index_t k, const T z) const
1211 {
1212  GISMO_ENSURE(m_data.mine().m_Acov_ori_mat.cols()!=0,"Is the metric initialized?");
1213  GISMO_ENSURE(m_data.mine().m_Acon_ori_mat.cols()!=0,"Is the metric initialized?");
1214  GISMO_ENSURE(m_data.mine().m_Bcov_ori_mat.cols()!=0,"Is the metric initialized?");
1215  GISMO_ENSURE(m_data.mine().m_ncov_ori_mat.cols()!=0,"Is the basis initialized?");
1216 
1217  GISMO_ENSURE(m_data.mine().m_acov_ori_mat.cols()!=0,"Is the basis initialized?");
1218  GISMO_ENSURE(m_data.mine().m_acon_ori_mat.cols()!=0,"Is the basis initialized?");
1219  GISMO_ENSURE(m_data.mine().m_normal_ori_mat.cols()!=0,"Is the basis initialized?");
1220 
1221  gsMatrix<T> Acov_ori, Acon_ori, Bcov_ori, ncov_ori, Gcov_ori(3,3), Gcon_ori(3,3),
1222  acov_ori, acon_ori, normal(3,1), gcov_ori(3,3), gcon_ori(3,3);
1223 
1224  // Get metric information
1225  Acov_ori = m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2);
1226  Acon_ori = m_data.mine().m_Acon_ori_mat.reshapeCol(k,2,2);
1227  Bcov_ori = m_data.mine().m_Bcov_ori_mat.reshapeCol(k,2,2);
1228  ncov_ori = m_data.mine().m_ncov_ori_mat.reshapeCol(k,3,2);
1229 
1230  // Compute full metric
1231  Gcov_ori.setZero();
1232  Gcov_ori.block(0,0,2,2)= Acov_ori - 2.0 * z * Bcov_ori + z*z * ncov_ori.transpose()*ncov_ori;
1233  Gcov_ori(2,2) = 1.0;
1234  Gcon_ori = Gcov_ori.inverse();
1235 
1236  // Assign members
1237  m_data.mine().m_Acov_ori = Acov_ori;
1238  m_data.mine().m_Acon_ori = Acon_ori;
1239  m_data.mine().m_Bcov_ori = Bcov_ori;
1240  m_data.mine().m_ncov_ori = ncov_ori;
1241  m_data.mine().m_Gcov_ori = Gcov_ori;
1242  m_data.mine().m_Gcon_ori = Gcon_ori;
1243 
1244  // Get basis vectors
1245  acov_ori = m_data.mine().m_acov_ori_mat.reshapeCol(k,3,2);
1246  acon_ori = m_data.mine().m_acon_ori_mat.reshapeCol(k,3,2);
1247  normal = m_data.mine().m_normal_ori_mat.reshapeCol(k,3,1);
1248 
1249  // Compute g_cov
1250  gcov_ori.setZero();
1251  gcov_ori.leftCols(2) = acov_ori + z * ncov_ori;
1252  gcov_ori.col(2) = normal;
1253 
1254  // Compute g_con
1255  for (index_t c = 0; c!=3; c++)
1256  gcon_ori.col(c) = Gcon_ori(c,0) * gcov_ori.col(0) + Gcon_ori(c,1) * gcov_ori.col(1) + Gcon_ori(c,2) * gcov_ori.col(2);
1257 
1258  // Assign members
1259  m_data.mine().m_acov_ori = acov_ori;
1260  m_data.mine().m_acon_ori = acov_ori;
1261  m_data.mine().m_gcov_ori = gcov_ori;
1262  m_data.mine().m_gcon_ori = gcon_ori;
1263 
1264  // // create a Linearised covariant tensor for SvK models
1265  // m_data.mine().m_Gcov_ori_L = m_data.mine().m_Gcov_ori;
1266  // m_data.mine().m_Gcov_ori.block(0,0,2,2) -= z*z * m_data.mine().m_ncov_ori.transpose()*m_data.mine().m_ncov_ori;
1267 }
1268 
1269 template <short_t dim, class T>
1270 template <short_t _dim>
1271 typename std::enable_if<_dim==2, void>::type
1272 gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed_impl(const index_t k, const T z) const
1273 {
1274  GISMO_ENSURE(m_data.mine().m_Acov_ori_mat.cols()!=0,"Is the metric initialized?");
1275  GISMO_ENSURE(m_data.mine().m_Acon_ori_mat.cols()!=0,"Is the metric initialized?");
1276  GISMO_ENSURE(m_data.mine().m_Bcov_ori_mat.cols()!=0,"Is the metric initialized?");
1277  GISMO_ENSURE(m_data.mine().m_ncov_ori_mat.cols()!=0,"Is the basis initialized?");
1278 
1279  GISMO_ENSURE(m_data.mine().m_acov_ori_mat.cols()!=0,"Is the basis initialized?");
1280  GISMO_ENSURE(m_data.mine().m_acon_ori_mat.cols()!=0,"Is the basis initialized?");
1281 
1282  gsMatrix<T> Acov_ori, Acon_ori, Gcov_ori(3,3), Gcon_ori(3,3),
1283  acov_ori, acon_ori, normal(3,1), gcov_ori(3,3), gcon_ori(3,3);
1284 
1285  // Get metric information
1286  Acov_ori = m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2);
1287  Acon_ori = m_data.mine().m_Acov_ori_mat.reshapeCol(k,2,2);
1288 
1289  // Compute full metric
1290  Gcov_ori.setZero();
1291  Gcov_ori.block(0,0,2,2)= Acov_ori;;
1292  Gcov_ori(2,2) = 1.0;
1293  Gcon_ori = Gcov_ori.inverse();
1294 
1295  // Assign members
1296  m_data.mine().m_Acov_ori = Acov_ori;
1297  m_data.mine().m_Acon_ori = Acon_ori;
1298  m_data.mine().m_Bcov_ori = m_data.mine().m_Bcov_ori_mat.reshapeCol(k,2,2);
1299  m_data.mine().m_ncov_ori = m_data.mine().m_ncov_ori_mat.reshapeCol(k,2,2);
1300  m_data.mine().m_Gcov_ori = Gcov_ori;
1301  m_data.mine().m_Gcon_ori = Gcon_ori;
1302 
1303  // Get basis vectors
1304  acov_ori = m_data.mine().m_acov_ori_mat.reshapeCol(k,2,2);
1305  acon_ori = m_data.mine().m_acon_ori_mat.reshapeCol(k,2,2);
1306  normal << 0,0,1;
1307 
1308  // Compute g_cov
1309  gcov_ori.setZero();
1310  gcov_ori.block(0,0,2,2) = acov_ori;
1311  gcov_ori.col(2) = normal;
1312 
1313  // Compute g_con
1314  for (index_t c = 0; c!=3; c++)
1315  gcon_ori.col(c) = Gcon_ori(c,0) * gcov_ori.col(0) + Gcon_ori(c,1) * gcov_ori.col(1) + Gcon_ori(c,2) * gcov_ori.col(2);
1316 
1317  // Assign members
1318  m_data.mine().m_acov_ori = acov_ori;
1319  m_data.mine().m_acon_ori = acon_ori;
1320  m_data.mine().m_gcov_ori = gcov_ori;
1321  m_data.mine().m_gcon_ori = gcon_ori;
1322 
1323  // // create a Linearised covariant tensor for SvK models
1324  // m_data.mine().m_Gcov_ori_L = m_data.mine().m_Gcov_ori;
1325  // m_data.mine().m_Gcov_ori.block(0,0,2,2) -= z*z * m_data.mine().m_ncov_ori.transpose()*m_data.mine().m_ncov_ori;
1326 }
1327 
1328 //--------------------------------------------------------------------------------------------------------------------------------------
1329 
1330 template <short_t dim, class T>
1331 std::pair<gsVector<T>,gsMatrix<T>> gsMaterialMatrixBaseDim<dim,T>::_evalStretch(const gsMatrix<T> & C, const gsMatrix<T> & gcon_ori) const
1332 {
1333  gsVector<T> stretches;
1334  gsMatrix<T> stretchvec;
1335  std::pair<gsVector<T>,gsMatrix<T>> result;
1336  stretches.resize(3,1); stretches.setZero();
1337  stretchvec.resize(3,3); stretchvec.setZero();
1338 
1339  typename gsEigen::SelfAdjointEigenSolver< typename gsMatrix<T>::Base > eigSolver;
1340 
1341  GISMO_ENSURE(m_data.mine().m_gcon_ori.cols()!=0,"Is the basis initialized?");
1342 
1343  gsMatrix<T> B(3,3);
1344  B.setZero();
1345  for (index_t k = 0; k != 2; k++)
1346  for (index_t l = 0; l != 2; l++)
1347  B += C(k,l) * gcon_ori.col(k) * gcon_ori.col(l).transpose();
1348 
1349  gsMatrix<T> eigVectors(3,3);
1350  gsMatrix<T> eigValues(3,1);
1351  if (math::abs(((B.block(0,0,2,2).diagonal()).array()-1).sum()) < 10*std::numeric_limits<T>::epsilon())
1352  {
1353  eigValues<<0,1,1; // eigenvalues when B = [1 0 0; 0 1 0; 0 0 0]
1354  eigVectors<<0,0,1,
1355  0,1,0,
1356  1,0,0;
1357  }
1358  else
1359  {
1360  eigSolver.compute(B);
1361  eigValues = eigSolver.eigenvalues();
1362  eigVectors= eigSolver.eigenvectors();
1363  }
1364 
1365 
1366  stretchvec.leftCols(2) = eigVectors.rightCols(2);
1367  stretchvec.col(2) = gcon_ori.col(2); // replace with: stretchvec.col(0).template head<3>().cross(stretchvec.col(1).template head<3>())
1368  stretches.block(0,0,2,1) = eigValues.block(1,0,2,1); // the eigenvalues are a 3x1 matrix, so we need to use matrix block-operations
1369 
1370  // m_data.mine().m_stretches.at(2) = 1/m_data.mine().m_J0_sq;
1371  stretches.at(2) = C(2,2);
1372 
1373  for (index_t k=0; k!=3; k++)
1374  stretches.at(k) = math::sqrt(stretches.at(k));
1375 
1376  result.first = stretches;
1377  result.second = stretchvec;
1378 
1379  return result;
1380 }
1381 
1382 
1383 template <short_t dim, class T>
1384 std::pair<gsVector<T>,gsMatrix<T>> gsMaterialMatrixBaseDim<dim,T>::_evalPStress(const gsMatrix<T> & S) const
1385 {
1386  gsVector<T> pstresses;
1387  gsMatrix<T> pstressvec;
1388  std::pair<gsVector<T>,gsMatrix<T>> result;
1389  pstresses.resize(2,1); pstresses.setZero();
1390  pstressvec.resize(3,2); pstressvec.setZero();
1391 
1392  typename gsEigen::SelfAdjointEigenSolver< typename gsMatrix<T>::Base > eigSolver;
1393 
1394  gsMatrix<T> B(3,3);
1395  B.setZero();
1396  for (index_t k = 0; k != 2; k++)
1397  for (index_t l = 0; l != 2; l++)
1398  B += S(k,l) * m_data.mine().m_gcov_ori.col(k) * m_data.mine().m_gcov_ori.col(l).transpose();
1399 
1400  eigSolver.compute(B);
1401 
1402  index_t zeroIdx = -1;
1403  T tol = 1e-14;
1404  T max = eigSolver.eigenvalues().array().abs().maxCoeff();
1405  max = (max==0) ? 1 : max;
1406  for (index_t k=0; k!=3; k++)
1407  zeroIdx = std::abs(eigSolver.eigenvalues()[k] ) / max < tol ? k : zeroIdx;
1408 
1409  GISMO_ASSERT(zeroIdx!=-1,"No zero found?");
1410 
1411  index_t count = 0;
1412 
1413  for (index_t k=0; k!=3; k++)
1414  {
1415  if (k==zeroIdx) continue;
1416  pstressvec.col(count) = eigSolver.eigenvectors().col(k);
1417  pstresses(count,0) = eigSolver.eigenvalues()(k,0);
1418  count++;
1419  }
1420 
1421  result.first = pstresses;
1422  result.second = pstressvec;
1423 
1424  return result;
1425 }
1426 
1427 template <short_t dim, class T>
1428 std::pair<gsVector<T>,gsMatrix<T>> gsMaterialMatrixBaseDim<dim,T>::_evalPStrain(const gsMatrix<T> & S) const
1429 {
1430  gsVector<T> pstrains;
1431  gsMatrix<T> pstrainvec;
1432  std::pair<gsVector<T>,gsMatrix<T>> result;
1433  pstrains.resize(3,1); pstrains.setZero();
1434  pstrainvec.resize(3,3); pstrainvec.setZero();
1435 
1436  typename gsEigen::SelfAdjointEigenSolver< typename gsMatrix<T>::Base > eigSolver;
1437 
1438  gsMatrix<T> B(3,3);
1439  B.setZero();
1440  for (index_t k = 0; k != 2; k++)
1441  for (index_t l = 0; l != 2; l++)
1442  B += S(k,l) * m_data.mine().m_gcon_ori.col(k) * m_data.mine().m_gcon_ori.col(l).transpose();
1443 
1444  eigSolver.compute(B);
1445 
1446  index_t zeroIdx = -1;
1447  T tol = 1e-14;
1448  T max = eigSolver.eigenvalues().array().abs().maxCoeff();
1449  max = (max==0) ? 1 : max;
1450  for (index_t k=0; k!=3; k++)
1451  zeroIdx = std::abs(eigSolver.eigenvalues()[k] ) / max < tol ? k : zeroIdx;
1452 
1453  GISMO_ASSERT(zeroIdx!=-1,"No zero found?");
1454 
1455  index_t count = 0;
1456  pstrainvec.col(2) = m_data.mine().m_gcon_ori.col(2);
1457  pstrains(2,0) = S(2,2);
1458 
1459  for (index_t k=0; k!=3; k++)
1460  {
1461  if (k==zeroIdx) continue;
1462  pstrainvec.col(count) = eigSolver.eigenvectors().col(k);
1463  pstrains(count,0) = eigSolver.eigenvalues()(k,0);
1464  count++;
1465  }
1466 
1467  result.first = pstrains;
1468  result.second = pstrainvec;
1469 
1470  return result;
1471 }
1472 
1473 //--------------------------------------------------------------------------------------------------------------------------------------
1474 
1475 template <short_t dim, class T>
1476 void gsMaterialMatrixBaseDim<dim,T>::_computeStretch(const gsMatrix<T> & C, const gsMatrix<T> & gcon_ori) const
1477 {
1478  std::pair<gsVector<T>,gsMatrix<T>> result = _evalStretch(C,gcon_ori);
1479  m_data.mine().m_stretches = result.first;
1480  m_data.mine().m_stretchvec = result.second;
1481 }
1482 
1483 template <short_t dim, class T>
1484 void gsMaterialMatrixBaseDim<dim,T>::_computePStress(const gsMatrix<T> & S) const
1485 {
1486  std::pair<gsVector<T>,gsMatrix<T>> result = _evalPStress(S);
1487  m_data.mine().m_pstress = result.first;
1488  m_data.mine().m_pstressvec = result.second;
1489 }
1490 
1491 template <short_t dim, class T>
1492 void gsMaterialMatrixBaseDim<dim,T>::_computePStrain(const gsMatrix<T> & E) const
1493 {
1494  std::pair<gsVector<T>,gsMatrix<T>> result = _evalPStrain(E);
1495  m_data.mine().m_pstrain = result.first;
1496  m_data.mine().m_pstrainvec = result.second;
1497 }
1498 
1499 //--------------------------------------------------------------------------------------------------------------------------------------
1500 
1501 template <short_t dim, class T>
1502 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_transformation(const gsMatrix<T> & basis1, const gsMatrix<T> & basis2) const
1503 {
1504  // Transformation of a quantity FROM basis2 TO basis1
1505  gsMatrix<T> Tmat(3,3);
1506 
1507  for (index_t i = 0; i!=2; i++)
1508  for (index_t j = 0; j!=2; j++)
1509  Tmat(i,j) = math::pow(basis2.col(i).dot(basis1.col(j)),2);
1510 
1511  Tmat(2,0) = basis2.col(1).dot(basis1.col(0)) * basis2.col(0).dot(basis1.col(0));
1512  Tmat(2,1) = basis2.col(0).dot(basis1.col(1)) * basis2.col(1).dot(basis1.col(1));
1513 
1514  Tmat(0,2) = 2*basis2.col(0).dot(basis1.col(1)) * basis2.col(0).dot(basis1.col(0));
1515  Tmat(1,2) = 2*basis2.col(1).dot(basis1.col(0)) * basis2.col(1).dot(basis1.col(1));
1516 
1517  Tmat(2,2) = basis2.col(0).dot(basis1.col(1)) * basis2.col(1).dot(basis1.col(0))
1518  +basis2.col(0).dot(basis1.col(0)) * basis2.col(1).dot(basis1.col(1));
1519 
1520  return Tmat;
1521 }
1522 
1523 } // end namespace
gismo::gsMaterialMatrixBaseDim::_getAcon_def
gsMatrix< T > _getAcon_def(index_t k, T z) const
Returns the contravariant a tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:658
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition: gsMaterialMatrixUtils.h:98
gsMaterialMatrixLinear.h
Provides linear material matrices.
gismo::gsMaterialMatrixBaseDim::_getgcov_ori_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getgcov_ori_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:969
gismo::gsMaterialMatrixBaseDim::_getBcov_ori
gsMatrix< T > _getBcov_ori(index_t k, T z) const
Returns the covariant b tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:739
gismo::gsMaterialMatrixBaseDim::eval3D_pstretchDir
virtual gsMatrix< T > eval3D_pstretchDir(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:222
gismo::gsMaterialMatrixBaseDim::defaultOptions
virtual void defaultOptions() override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:38
gismo::gsMaterialMatrixBaseDim::_computeMetricDeformed
void _computeMetricDeformed(const index_t patch, const gsMatrix< T > &u) const
Computes metric quantities on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:417
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
gismo::gsMaterialMatrixBaseDim::_evalStretch
std::pair< gsVector< T >, gsMatrix< T > > _evalStretch(const gsMatrix< T > &C, const gsMatrix< T > &gcon_ori) const
Computes the stretch given deformation tensor C, into a pair.
Definition: gsMaterialMatrixBaseDim.hpp:1331
gismo::gsMaterialMatrixBaseDim::_getgcon_def
gsMatrix< T > _getgcon_def(index_t k, T z) const
Returns the contravariant basis vector g on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:710
gismo::NEED_DERIV
Gradient of the object.
Definition: gsForwardDeclarations.h:73
gismo::gsMaterialMatrixBaseDim::_getGcov_def_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getGcov_def_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:837
gismo::gsMaterialMatrixBaseDim::_getBcov_def
gsMatrix< T > _getBcov_def(index_t k, T z) const
Returns the covariant b tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:665
gismo::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition: gsLinearAlgebra.h:126
gismo::gsMaterialMatrixBaseDim::eval3D_pstretch
virtual gsMatrix< T > eval3D_pstretch(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:200
gismo::gsMaterialMatrixBaseDim::_getAcon_ori
gsMatrix< T > _getAcon_ori(index_t k, T z) const
Returns the contravariant a tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:732
gismo::gsMaterialMatrixBaseDim::eval3D_pstrain
virtual gsMatrix< T > eval3D_pstrain(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:244
gismo::gsMaterialMatrixBaseDim::eval3D_strain
virtual gsMatrix< T > eval3D_strain(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:271
gismo::gsMaterialMatrixBaseDim::_getacon_def
gsMatrix< T > _getacon_def(index_t k, T z) const
Returns the contravariant basis vector a on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:697
gismo::gsMaterialMatrixBaseDim::_getncov_def
gsMatrix< T > _getncov_def(index_t k, T z) const
Returns the covariant n tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:671
gismo::gsMaterialMatrixBaseDim::eval3D_con2cart
virtual gsMatrix< T > eval3D_con2cart(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:137
gismo::gsMaterialMatrixBaseDim::_getGcon_def
gsMatrix< T > _getGcon_def(index_t k, T z) const
Returns the contravariant metric tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:683
gismo::gsMaterialMatrixBaseDim::eval3D_spec2cov
virtual gsMatrix< T > eval3D_spec2cov(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:72
gismo::gsMaterialMatrixBaseDim::_getacon_ori
gsMatrix< T > _getacon_ori(index_t k, T z) const
Returns the contravariant basis vector a on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:771
gismo::gsMaterialMatrixBaseDim::_evalPStress
std::pair< gsVector< T >, gsMatrix< T > > _evalPStress(const gsMatrix< T > &S) const
Computes the principal stress given stress tensor S, into a pair.
Definition: gsMaterialMatrixBaseDim.hpp:1384
gismo::gsMaterialMatrixBaseDim::_getGcov_ori
gsMatrix< T > _getGcov_ori(index_t k, T z) const
Returns the covariant metric tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:751
gismo::gsMapData
the gsMapData is a cache of pre-computed function (map) values.
Definition: gsFuncData.h:324
gismo::gsMaterialMatrixBaseDim::_getBcov_def_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getBcov_def_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:803
gismo::gsMaterialMatrixBaseDim::density_into
virtual void density_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:44
gismo::gsMaterialMatrixBaseDim::eval3D_cov2cart
virtual gsMatrix< T > eval3D_cov2cart(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:116
gsGeometry.h
Provides declaration of Geometry abstract interface.
gismo::NEED_DERIV2
Second derivatives.
Definition: gsForwardDeclarations.h:80
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition: gsFunction.h:59
gismo::gsMaterialMatrixBaseDim::_computeMetricUndeformed_impl
std::enable_if< _dim==2, void >::type _computeMetricUndeformed_impl(const index_t patch, const gsMatrix< T > &u) const
Implementation of _getMetric for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:606
gismo::gsMaterialMatrixBaseDim::parameters_into
virtual void parameters_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:157
gismo::gsMaterialMatrixBaseDim::_computePStrain
void _computePStrain(const gsMatrix< T > &C) const
Computes the stretch given deformation tensor C, into class members m_stretches and m_stretchDirs...
Definition: gsMaterialMatrixBaseDim.hpp:1492
gismo::gsMaterialMatrixBaseDim::_getacov_def
gsMatrix< T > _getacov_def(index_t k, T z) const
Returns the covariant basis vector a on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:690
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gismo::gsMaterialMatrixBaseDim::_getgcov_def
gsMatrix< T > _getgcov_def(index_t k, T z) const
Returns the covariant basis vector g on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:704
gsMaterialMatrixBaseDim.h
Base class with dimension in template; used for metric computations.
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::gsMaterialMatrixBaseDim::_evalPStrain
std::pair< gsVector< T >, gsMatrix< T > > _evalPStrain(const gsMatrix< T > &C) const
Computes the principal strain given deformation tensor C, into a pair.
Definition: gsMaterialMatrixBaseDim.hpp:1428
gismo::gsMaterialMatrixBaseDim::thickness_into
virtual void thickness_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:61
gismo::gsMaterialMatrixBaseDim::_getMetricDeformed
void _getMetricDeformed(const index_t k, const T z) const
Gets metric quantities on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:1073
gismo::gsMaterialMatrixBaseDim::_computeMetricDeformed_impl
std::enable_if< _dim==2, void >::type _computeMetricDeformed_impl(const index_t patch, const gsMatrix< T > &u) const
Implementation of _computeMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:488
gismo::gsMaterialMatrixBaseDim::eval3D_tensionfield
virtual gsMatrix< T > eval3D_tensionfield(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:296
gsWarn
#define gsWarn
Definition: gsDebug.h:50
gismo::gsMaterialMatrixBaseDim::_getGcon_ori
gsMatrix< T > _getGcon_ori(index_t k, T z) const
Returns the contravariant metric tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:757
gismo::NEED_NORMAL
Normal vector of the object.
Definition: gsForwardDeclarations.h:85
gismo::gsMaterialMatrixBaseDim::_getMetricUndeformed
void _getMetricUndeformed(const index_t k, const T z) const
Gets metric quantities on the undeformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:1202
gismo::gsMaterialMatrixBaseDim::_computeMetricUndeformed
void _computeMetricUndeformed(const index_t patch, const gsMatrix< T > &u) const
Computes metric quantities on the undeformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:535
gismo::gsMaterialMatrixBaseDim::_getMetricDeformed_impl
std::enable_if< _dim==2, void >::type _getMetricDeformed_impl(const index_t k, const T z) const
Implementation of _getMetricDeformed for surface geometries.
Definition: gsMaterialMatrixBaseDim.hpp:1141
gismo::gsMaterialMatrixBaseDim::_getBcov_ori_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getBcov_ori_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:901
gismo::gsMaterialMatrixBaseDim::_getncov_def_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getncov_def_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:820
gismo::gsMaterialMatrixBaseDim::_getgcov_def_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getgcov_def_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:871
gismo::gsMaterialMatrixBaseDim::_getAcov_ori
gsMatrix< T > _getAcov_ori(index_t k, T z) const
Returns the covariant a tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:725
gismo::gsMaterialMatrixBaseDim::_computePStress
void _computePStress(const gsMatrix< T > &C) const
Computes the principal stresses of a given stress tensor S, into class members m_pstress and m_pstres...
Definition: gsMaterialMatrixBaseDim.hpp:1484
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::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition: gsVector.h:177
gsMaterialMatrixNonlinear.h
Provides hyperelastic material matrices.
gismo::gsMaterialMatrixBaseDim
This class defines the base class for material matrices.
Definition: gsMaterialMatrixBaseDim.h:35
gismo::gsMaterialMatrixBaseDim::_transformation
gsMatrix< T > _transformation(const gsMatrix< T > &basis1, const gsMatrix< T > &basis2) const
Computes the stretch given deformation tensor C, into class members m_stretches and m_stretchDirs...
Definition: gsMaterialMatrixBaseDim.hpp:1502
gismo::gsMaterialMatrixBaseDim::_computeStretch
void _computeStretch(const gsMatrix< T > &C, const gsMatrix< T > &gcon_ori) const
Computes the stretch given deformation tensor C, into class members m_stretches and m_stretchDirs...
Definition: gsMaterialMatrixBaseDim.hpp:1476
gismo::NEED_JACOBIAN
Jacobian of the object.
Definition: gsForwardDeclarations.h:75
gismo::gsMaterialMatrixBaseDim::_getgcov_ori
gsMatrix< T > _getgcov_ori(index_t k, T z) const
Returns the covariant basis vector g on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:778
gismo::gsMaterialMatrixBaseDim::_getGcov_def
gsMatrix< T > _getGcov_def(index_t k, T z) const
Returns the covariant metric tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:677
gismo::gsMaterialMatrixBaseDim::eval3D_spec2con
virtual gsMatrix< T > eval3D_spec2con(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:94
gismo::gsMaterialMatrixBaseDim::_getncov_ori_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getncov_ori_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:918
gismo::NEED_VALUE
Value of the object.
Definition: gsForwardDeclarations.h:72
gismo::gsMaterialMatrixBaseDim::_getacov_ori
gsMatrix< T > _getacov_ori(index_t k, T z) const
Returns the covariant basis vector a on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:764
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
gismo::gsMaterialMatrixBaseDim::_getMetricUndeformed_impl
std::enable_if< _dim==2, void >::type _getMetricUndeformed_impl(const index_t k, const T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:1272
gismo::gsMapData::points
gsMatrix< T > points
input (parametric) points
Definition: gsFuncData.h:348
gismo::gsMaterialMatrixBaseDim::_getAcov_def
gsMatrix< T > _getAcov_def(index_t k, T z) const
Returns the covariant a tensor on the deformed geometry.
Definition: gsMaterialMatrixBaseDim.hpp:651
gismo::gsMaterialMatrixBaseDim::_getGcov_ori_impl
std::enable_if< _dim==2, gsMatrix< T > >::type _getGcov_ori_impl(index_t k, T z) const
Implementation of _getMetricUndeformed for planar geometries.
Definition: gsMaterialMatrixBaseDim.hpp:935
gismo::gsMaterialMatrixBaseDim::eval3D_deformation
virtual gsMatrix< T > eval3D_deformation(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixBaseDim.hpp:176
gismo::gsMaterialMatrixBaseDim::_getgcon_ori
gsMatrix< T > _getgcon_ori(index_t k, T z) const
Returns the contravariant basis vector g on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:784
gismo::gsMaterialMatrixBaseDim::_getncov_ori
gsMatrix< T > _getncov_ori(index_t k, T z) const
Returns the covariant n tensor on the original geometry.
Definition: gsMaterialMatrixBaseDim.hpp:745
gismo::gsMaterialMatrixBaseDim::_getMetric
void _getMetric(const index_t k, const T z) const
Gets metric quantities on the deformed and undeformed geometries.
Definition: gsMaterialMatrixBaseDim.hpp:1021