1 /** @file gsMaterialMatrix.hpp
2
3 @brief Provides hyperelastic material matrices
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s):
12 H.M. Verhelst (2019-..., TU Delft)
13 A. Mantzaflaris (2019-..., Inria)
14 */
15
16 /*
17 To Do [updated 16-06-2020]:
18 - Make beta (compressible materials) and material parameters universal for all integration points over the thickness. So get them out of the _dPsi functions etc and move them into the integration loops as global variables.
19
20 */
21
22
23
24 #pragma once
25
26 #include <gsKLShell/gsMaterialMatrixBaseDim.h>
27
28 namespace gismo
29 {
30
31 //--------------------------------------------------------------------------------------------------------------------------------------
32
33 template <short_t dim, class T>
_computeMetricDeformed(const index_t patch,const gsMatrix<T> & u,bool basis) const34 void gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed(const index_t patch, const gsMatrix<T> & u, bool basis) const
35 {
36 m_map_def.mine().flags = m_map.mine().flags;
37 m_map_def.mine().points = u;
38 static_cast<const gsFunction<T>&>(Base::m_defpatches->piece(patch)).computeMap(m_map_def); // the piece(0) here implies that if you call class.eval_into, it will be evaluated on piece(0). Hence, call class.piece(k).eval_into()
39
40 _computeMetricDeformed_impl<dim>(patch,u,basis);
41 }
42
43 template <short_t dim, class T>
44 template <short_t _dim>
45 typename std::enable_if<_dim==3, void>::type
_computeMetricDeformed_impl(const index_t patch,const gsMatrix<T> & u,bool basis) const46 gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed_impl(const index_t patch, const gsMatrix<T> & u, bool basis) const
47 {
48 gsMatrix<T> deriv2;
49 gsMatrix<T,3,1> normal;
50 gsMatrix<T,2,2> mixedB;
51 gsMatrix<T,2,2> tmp;
52
53 m_Acov_def_mat.resize(4,m_map_def.mine().points.cols()); m_Acov_def_mat.setZero();
54 m_Acon_def_mat.resize(4,m_map_def.mine().points.cols()); m_Acon_def_mat.setZero();
55 m_Bcov_def_mat.resize(4,m_map_def.mine().points.cols()); m_Bcov_def_mat.setZero();
56
57 m_acov_def_mat.resize(2*3,m_map_def.mine().points.cols()); m_acov_def_mat.setZero();
58 if (basis)
59 {
60 m_acon_def_mat.resize(2*3,m_map_def.mine().points.cols()); m_acon_def_mat.setZero();
61 m_ncov_def_mat.resize(2*3,m_map_def.mine().points.cols()); m_ncov_def_mat.setZero();
62 }
63
64 for (index_t k=0; k!= m_map_def.mine().points.cols(); k++)
65 {
66 m_acov_def_mat.reshapeCol(k,3,2) = m_map_def.mine().jacobian(k);
67 gsAsMatrix<T,Dynamic,Dynamic> acov = m_acov_def_mat.reshapeCol(k,3,2);
68 acov = m_map_def.mine().jacobian(k);
69
70 tmp = acov.transpose() * acov;
71 m_Acov_def_mat.reshapeCol(k,2,2) = tmp;
72 m_Acon_def_mat.reshapeCol(k,2,2) = tmp.inverse();
73
74 gsAsMatrix<T,Dynamic,Dynamic> metricAcon = m_Acon_def_mat.reshapeCol(k,2,2);
75
76 // Construct metric tensor b = [d11c*n, d12c*n ; d21c*n, d22c*n]
77 deriv2 = m_map_def.mine().deriv2(k);
78 deriv2.resize(3,3);
79 normal = m_map_def.mine().normal(k).normalized();
80
81 tmp(0,0) = deriv2.row(0).dot(normal);
82 tmp(1,1) = deriv2.row(1).dot(normal);
83 tmp(0,1) = tmp(1,0) = deriv2.row(2).dot(normal);
84
85 m_Bcov_def_mat.reshapeCol(k,2,2) = tmp;
86 gsAsMatrix<T,Dynamic,Dynamic> metricBcov = m_Bcov_def_mat.reshapeCol(k,2,2);
87
88 // Construct basis
89 if (basis)
90 {
91 gsAsMatrix<T,Dynamic,Dynamic> acon = m_acon_def_mat.reshapeCol(k,3,2);
92 for (index_t i=0; i < 2; i++)
93 acon.col(i) = metricAcon(i,0)*acov.col(0) + metricAcon(i,1)*acov.col(1);
94
95 // Mixed tensor
96 for (index_t i=0; i < 2; i++)
97 for (index_t j=0; j < 2; j++)
98 mixedB(i,j) = metricAcon(i,0)*metricBcov(0,j) + metricAcon(i,1)*metricBcov(1,j);
99
100 gsAsMatrix<T,Dynamic,Dynamic> ncov = m_ncov_def_mat.reshapeCol(k,3,2);
101 for (index_t i=0; i < 2; i++)
102 ncov.col(i) = -mixedB(0,i)*acov.col(0) -mixedB(1,i)*acov.col(1);
103 }
104 }
105 }
106
107 template <short_t dim, class T>
108 template <short_t _dim>
109 typename std::enable_if<_dim==2, void>::type
_computeMetricDeformed_impl(const index_t patch,const gsMatrix<T> & u,bool basis) const110 gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed_impl(const index_t patch, const gsMatrix<T> & u, bool basis) const
111 {
112 gsMatrix<T,2,2> tmp;
113
114 m_Acov_def_mat.resize(4,m_map_def.mine().points.cols()); m_Acov_def_mat.setZero();
115 m_Acon_def_mat.resize(4,m_map_def.mine().points.cols()); m_Acon_def_mat.setZero();
116
117 m_acov_def_mat.resize(2*2,m_map_def.mine().points.cols()); m_acov_def_mat.setZero();
118 if (basis)
119 {
120 m_acon_def_mat.resize(2*2,m_map_def.mine().points.cols()); m_acon_def_mat.setZero();
121 }
122
123 for (index_t k=0; k!= m_map_def.mine().points.cols(); k++)
124 {
125 gsAsMatrix<T,Dynamic,Dynamic> acov = m_acov_def_mat.reshapeCol(k,2,2);
126 acov = m_map_def.mine().jacobian(k);
127
128 tmp = acov.transpose() * acov;
129 m_Acov_def_mat.reshapeCol(k,2,2) = tmp;
130 m_Acon_def_mat.reshapeCol(k,2,2) = tmp.inverse();
131
132 gsAsMatrix<T,Dynamic,Dynamic> metricAcon = m_Acon_def_mat.reshapeCol(k,2,2);
133
134 // Construct basis
135 if (basis)
136 {
137 gsAsMatrix<T,Dynamic,Dynamic> acon = m_acon_def_mat.reshapeCol(k,2,2);
138 for (index_t i=0; i < 2; i++)
139 acon.col(i) = metricAcon(i,0)*acov.col(0) + metricAcon(i,1)*acov.col(1);
140 }
141 }
142 }
143
144 //--------------------------------------------------------------------------------------------------------------------------------------
145
146 template <short_t dim, class T>
_computeMetricUndeformed(const index_t patch,const gsMatrix<T> & u,bool basis) const147 void gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed(const index_t patch, const gsMatrix<T> & u, bool basis) const
148 {
149 m_map.mine().flags = NEED_JACOBIAN | NEED_DERIV | NEED_NORMAL | NEED_VALUE | NEED_DERIV2;
150 m_map.mine().points = u;
151 static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(m_map); // the piece(0) here implies that if you call class.eval_into, it will be evaluated on piece(0). Hence, call class.piece(k).eval_into()
152
153 _computeMetricUndeformed_impl<dim>(patch,u,basis);
154 }
155
156 template <short_t dim, class T>
157 template <short_t _dim>
158 typename std::enable_if<_dim==3, void>::type
_computeMetricUndeformed_impl(const index_t patch,const gsMatrix<T> & u,bool basis) const159 gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed_impl(const index_t patch, const gsMatrix<T> & u, bool basis) const
160 {
161 gsMatrix<T> deriv2;
162 gsMatrix<T,3,1> normal;
163 gsMatrix<T,2,2> mixedB;
164 gsMatrix<T,2,2> tmp;
165
166 m_Acov_ori_mat.resize(4,m_map.mine().points.cols()); m_Acov_ori_mat.setZero();
167 m_Acon_ori_mat.resize(4,m_map.mine().points.cols()); m_Acon_ori_mat.setZero();
168 m_Bcov_ori_mat.resize(4,m_map.mine().points.cols()); m_Bcov_ori_mat.setZero();
169
170 m_acov_ori_mat.resize(2*3,m_map.mine().points.cols()); m_acov_ori_mat.setZero();
171 if (basis)
172 {
173 m_acon_ori_mat.resize(2*3,m_map.mine().points.cols()); m_acon_ori_mat.setZero();
174 m_ncov_ori_mat.resize(2*3,m_map.mine().points.cols()); m_ncov_ori_mat.setZero();
175 }
176
177 for (index_t k=0; k!= m_map.mine().points.cols(); k++)
178 {
179 m_acov_ori_mat.reshapeCol(k,3,2) = m_map.mine().jacobian(k);
180 gsAsMatrix<T,Dynamic,Dynamic> acov = m_acov_ori_mat.reshapeCol(k,3,2);
181 acov = m_map.mine().jacobian(k);
182
183 tmp = acov.transpose() * acov;
184 m_Acov_ori_mat.reshapeCol(k,2,2) = tmp;
185 m_Acon_ori_mat.reshapeCol(k,2,2) = tmp.inverse();
186
187 gsAsMatrix<T,Dynamic,Dynamic> metricAcon = m_Acon_ori_mat.reshapeCol(k,2,2);
188
189 // Construct metric tensor b = [d11c*n, d12c*n ; d21c*n, d22c*n]
190 deriv2 = m_map.mine().deriv2(k);
191 deriv2.resize(3,3);
192 normal = m_map.mine().normal(k).normalized();
193
194 tmp(0,0) = deriv2.row(0).dot(normal);
195 tmp(1,1) = deriv2.row(1).dot(normal);
196 tmp(0,1) = tmp(1,0) = deriv2.row(2).dot(normal);
197
198 m_Bcov_ori_mat.reshapeCol(k,2,2) = tmp;
199 gsAsMatrix<T,Dynamic,Dynamic> metricBcov = m_Bcov_ori_mat.reshapeCol(k,2,2);
200
201 // Construct basis
202 if (basis)
203 {
204 gsAsMatrix<T,Dynamic,Dynamic> acon = m_acon_ori_mat.reshapeCol(k,3,2);
205 for (index_t i=0; i < 2; i++)
206 acon.col(i) = metricAcon(i,0)*acov.col(0) + metricAcon(i,1)*acov.col(1);
207
208 // Mixed tensor
209 for (index_t i=0; i < 2; i++)
210 for (index_t j=0; j < 2; j++)
211 mixedB(i,j) = metricAcon(i,0)*metricBcov(0,j) + metricAcon(i,1)*metricBcov(1,j);
212
213 gsAsMatrix<T,Dynamic,Dynamic> ncov = m_ncov_ori_mat.reshapeCol(k,3,2);
214 for (index_t i=0; i < 2; i++)
215 ncov.col(i) = -mixedB(0,i)*acov.col(0) -mixedB(1,i)*acov.col(1);
216 }
217 }
218 }
219
220 template <short_t dim, class T>
221 template <short_t _dim>
222 typename std::enable_if<_dim==2, void>::type
_computeMetricUndeformed_impl(const index_t patch,const gsMatrix<T> & u,bool basis) const223 gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed_impl(const index_t patch, const gsMatrix<T> & u, bool basis) const
224 {
225 gsMatrix<T,2,2> tmp;
226
227 m_Acov_ori_mat.resize(4,m_map.mine().points.cols()); m_Acov_ori_mat.setZero();
228 m_Acon_ori_mat.resize(4,m_map.mine().points.cols()); m_Acon_ori_mat.setZero();
229
230 m_acov_ori_mat.resize(2*2,m_map.mine().points.cols()); m_acov_ori_mat.setZero();
231 if (basis)
232 {
233 m_acon_ori_mat.resize(2*2,m_map.mine().points.cols()); m_acon_ori_mat.setZero();
234 }
235
236 for (index_t k=0; k!= m_map.mine().points.cols(); k++)
237 {
238 gsAsMatrix<T,Dynamic,Dynamic> acov = m_acov_ori_mat.reshapeCol(k,2,2);
239 acov = m_map.mine().jacobian(k);
240
241 tmp = acov.transpose() * acov;
242 m_Acov_ori_mat.reshapeCol(k,2,2) = tmp;
243 m_Acon_ori_mat.reshapeCol(k,2,2) = tmp.inverse();
244
245 gsAsMatrix<T,Dynamic,Dynamic> metricAcon = m_Acon_ori_mat.reshapeCol(k,2,2);
246
247 // Construct basis
248 if (basis)
249 {
250 gsAsMatrix<T,Dynamic,Dynamic> acon = m_acon_ori_mat.reshapeCol(k,2,2);
251 for (index_t i=0; i < 2; i++)
252 acon.col(i) = metricAcon(i,0)*acov.col(0) + metricAcon(i,1)*acov.col(1);
253 }
254 }
255 }
256
257 //--------------------------------------------------------------------------------------------------------------------------------------
258
259 template <short_t dim, class T>
_getMetric(index_t k,T z,bool basis) const260 void gsMaterialMatrixBaseDim<dim,T>::_getMetric(index_t k, T z, bool basis) const
261 {
262 this->_getMetricDeformed(k,z,basis);
263 this->_getMetricUndeformed(k,z,basis);
264
265 T ratio = m_Gcov_def.determinant() / m_Gcov_ori.determinant();
266 GISMO_ENSURE(ratio > 0, "Jacobian determinant is negative! det(Gcov_def) = "<<m_Gcov_def.determinant()<<"; det(Gcov_ori) = "<<m_Gcov_ori.determinant());
267 m_J0_sq = ratio;
268 }
269
270 //--------------------------------------------------------------------------------------------------------------------------------------
271
272 template <short_t dim, class T>
_getMetricDeformed(index_t k,T z,bool basis) const273 void gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed(index_t k, T z, bool basis) const
274 {
275 _getMetricDeformed_impl<dim>(k,z,basis);
276 }
277
278 template <short_t dim, class T>
279 template <short_t _dim>
280 typename std::enable_if<_dim==3, void>::type
_getMetricDeformed_impl(index_t k,T z,bool basis) const281 gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed_impl(index_t k, T z, bool basis) const
282 {
283 GISMO_ENSURE(m_Acov_def_mat.cols()!=0,"Is the metric initialized?");
284 GISMO_ENSURE(m_Acon_def_mat.cols()!=0,"Is the metric initialized?");
285 GISMO_ENSURE(m_Bcov_def_mat.cols()!=0,"Is the metric initialized?");
286 if (basis)
287 {
288 GISMO_ENSURE(m_acov_def_mat.cols()!=0,"Is the basis initialized?");
289 GISMO_ENSURE(m_acon_def_mat.cols()!=0,"Is the basis initialized?");
290 GISMO_ENSURE(m_ncov_def_mat.cols()!=0,"Is the basis initialized?");
291 }
292
293 // metrics
294 m_Acov_def = m_Acov_def_mat.reshapeCol(k,2,2);
295 m_Acon_def = m_Acon_def_mat.reshapeCol(k,2,2);
296 m_Bcov_def = m_Bcov_def_mat.reshapeCol(k,2,2);
297 if (basis)
298 m_ncov_def = m_ncov_def_mat.reshapeCol(k,3,2);
299
300 // Compute full metric
301 m_Gcov_def.setZero();
302 m_Gcov_def.block(0,0,2,2)= m_Acov_def - 2.0 * z * m_Bcov_def + z*z * m_ncov_def.transpose()*m_ncov_def;
303 m_Gcov_def(2,2) = 1.0;
304 m_Gcon_def = m_Gcov_def.inverse();
305
306 if (!basis) return;
307 // Compute full basis
308 // basis vectors
309 m_acov_def = m_acov_def_mat.reshapeCol(k,3,2);
310 m_acon_def = m_acon_def_mat.reshapeCol(k,3,2);
311 // g
312 gsMatrix<T,3,1> normal = m_map_def.mine().normal(k).normalized();
313 m_gcov_def.leftCols(2) = m_acov_def + z * m_ncov_def;
314 m_gcov_def.col(2) = normal;
315
316 for (index_t c = 0; c!=3; c++)
317 {
318 m_gcon_def.col(c) = m_Gcon_def(c,0) * m_gcov_def.col(0)
319 + m_Gcon_def(c,1) * m_gcov_def.col(1)
320 + m_Gcon_def(c,2) * m_gcov_def.col(2);
321 }
322
323 // // create a Linearised covariant tensor for SvK models
324 // m_Gcov_def_L = m_Gcov_def;
325 // m_Gcov_def.block(0,0,2,2) -= z*z * m_ncov_def.transpose()*m_ncov_def;
326 }
327
328 template <short_t dim, class T>
329 template <short_t _dim>
330 typename std::enable_if<_dim==2, void>::type
_getMetricDeformed_impl(index_t k,T z,bool basis) const331 gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed_impl(index_t k, T z, bool basis) const
332 {
333 GISMO_ENSURE(m_Acov_def_mat.cols()!=0,"Is the metric initialized?");
334 GISMO_ENSURE(m_Acon_def_mat.cols()!=0,"Is the metric initialized?");
335
336 if (basis)
337 {
338 GISMO_ENSURE(m_acov_def_mat.cols()!=0,"Is the basis initialized?");
339 GISMO_ENSURE(m_acon_def_mat.cols()!=0,"Is the basis initialized?");
340 }
341
342 // metrics
343 m_Acov_def = m_Acov_def_mat.reshapeCol(k,2,2);
344 m_Acon_def = m_Acon_def_mat.reshapeCol(k,2,2);
345 // Compute full metric
346 m_Gcov_def.setZero();
347 m_Gcov_def.block(0,0,2,2)= m_Acov_def;
348 m_Gcov_def(2,2) = 1.0;
349 m_Gcon_def = m_Gcov_def.inverse();
350
351 if (!basis) return;
352 // Compute full basis
353 // basis vectors
354 m_acov_def = m_acov_def_mat.reshapeCol(k,2,2);
355 m_acon_def = m_acon_def_mat.reshapeCol(k,2,2);
356 // g
357 gsMatrix<T,3,1> normal;
358 normal << 0,0,1;
359 m_gcov_def.setZero();
360 m_gcov_def.block(0,0,2,2) = m_acov_def;
361 m_gcov_def.col(2) = normal;
362
363 for (index_t c = 0; c!=3; c++)
364 {
365 m_gcon_def.col(c) = m_Gcon_def(c,0) * m_gcov_def.col(0)
366 + m_Gcon_def(c,1) * m_gcov_def.col(1)
367 + m_Gcon_def(c,2) * m_gcov_def.col(2);
368 }
369
370 // // create a Linearised covariant tensor for SvK models
371 // m_Gcov_def_L = m_Gcov_def;
372 // m_Gcov_def.block(0,0,2,2) -= z*z * m_ncov_def.transpose()*m_ncov_def;
373 }
374
375 //--------------------------------------------------------------------------------------------------------------------------------------
376
377 template <short_t dim, class T>
_getMetricUndeformed(index_t k,T z,bool basis) const378 void gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed(index_t k, T z, bool basis) const
379 {
380 _getMetricUndeformed_impl<dim>(k,z,basis);
381 }
382
383 template <short_t dim, class T>
384 template <short_t _dim>
385 typename std::enable_if<_dim==3, void>::type
_getMetricUndeformed_impl(index_t k,T z,bool basis) const386 gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed_impl(index_t k, T z, bool basis) const
387 {
388 GISMO_ENSURE(m_Acov_ori_mat.cols()!=0,"Is the metric initialized?");
389 GISMO_ENSURE(m_Acon_ori_mat.cols()!=0,"Is the metric initialized?");
390 GISMO_ENSURE(m_Bcov_ori_mat.cols()!=0,"Is the metric initialized?");
391
392 if (basis)
393 {
394 GISMO_ENSURE(m_acov_ori_mat.cols()!=0,"Is the basis initialized?");
395 GISMO_ENSURE(m_acon_ori_mat.cols()!=0,"Is the basis initialized?");
396 GISMO_ENSURE(m_ncov_ori_mat.cols()!=0,"Is the basis initialized?");
397 }
398
399 // metrics
400 m_Acov_ori = m_Acov_ori_mat.reshapeCol(k,2,2);
401 m_Acon_ori = m_Acon_ori_mat.reshapeCol(k,2,2);
402 m_Bcov_ori = m_Bcov_ori_mat.reshapeCol(k,2,2);
403 if (basis)
404 m_ncov_ori = m_ncov_ori_mat.reshapeCol(k,3,2);
405
406 // Compute full metric
407 m_Gcov_ori.setZero();
408 m_Gcov_ori.block(0,0,2,2)= m_Acov_ori - 2.0 * z * m_Bcov_ori + z*z * m_ncov_ori.transpose()*m_ncov_ori;
409 m_Gcov_ori(2,2) = 1.0;
410 m_Gcon_ori = m_Gcov_ori.inverse();
411
412 if (!basis) return;
413 // Compute full basis
414 // basis vectors
415 m_acov_ori = m_acov_ori_mat.reshapeCol(k,3,2);
416 m_acon_ori = m_acon_ori_mat.reshapeCol(k,3,2);
417 // g
418 gsMatrix<T,3,1> normal = m_map.mine().normal(k).normalized();
419 m_gcov_ori.block(0,0,3,2) = m_acov_ori + z * m_ncov_ori;
420 m_gcov_ori.col(2) = normal;
421 for (index_t c = 0; c!=3; c++)
422 {
423 m_gcon_ori.col(c) = m_Gcon_ori(c,0) * m_gcov_ori.col(0)
424 + m_Gcon_ori(c,1) * m_gcov_ori.col(1)
425 + m_Gcon_ori(c,2) * m_gcov_ori.col(2);
426 }
427
428 // // create a Linearised covariant tensor for SvK models
429 // m_Gcov_ori_L = m_Gcov_ori;
430 // m_Gcov_ori.block(0,0,2,2) -= z*z * m_ncov_ori.transpose()*m_ncov_ori;
431 }
432
433 template <short_t dim, class T>
434 template <short_t _dim>
435 typename std::enable_if<_dim==2, void>::type
_getMetricUndeformed_impl(index_t k,T z,bool basis) const436 gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed_impl(index_t k, T z, bool basis) const
437 {
438 GISMO_ENSURE(m_Acov_ori_mat.cols()!=0,"Is the metric initialized?");
439 GISMO_ENSURE(m_Acon_ori_mat.cols()!=0,"Is the metric initialized?");
440
441 if (basis)
442 {
443 GISMO_ENSURE(m_acov_ori_mat.cols()!=0,"Is the basis initialized?");
444 GISMO_ENSURE(m_acon_ori_mat.cols()!=0,"Is the basis initialized?");
445 }
446
447 // metrics
448 m_Acov_ori = m_Acov_ori_mat.reshapeCol(k,2,2);
449 m_Acon_ori = m_Acon_ori_mat.reshapeCol(k,2,2);
450 // Compute full metric
451 m_Gcov_ori.setZero();
452 m_Gcov_ori.block(0,0,2,2)= m_Acov_ori;
453 m_Gcov_ori(2,2) = 1.0;
454 m_Gcon_ori = m_Gcov_ori.inverse();
455
456 if (!basis) return;
457 // Compute full basis
458 // basis vectors
459 m_acov_ori = m_acov_ori_mat.reshapeCol(k,2,2);
460 m_acon_ori = m_acon_ori_mat.reshapeCol(k,2,2);
461 // g
462 gsMatrix<T,3,1> normal;
463 normal << 0,0,1;
464 m_gcov_ori.setZero();
465 m_gcov_ori.block(0,0,2,2) = m_acov_ori;
466 m_gcov_ori.col(2) = normal;
467
468 for (index_t c = 0; c!=3; c++)
469 {
470 m_gcon_ori.col(c) = m_Gcon_ori(c,0) * m_gcov_ori.col(0)
471 + m_Gcon_ori(c,1) * m_gcov_ori.col(1)
472 + m_Gcon_ori(c,2) * m_gcov_ori.col(2);
473 }
474
475 // // create a Linearised covariant tensor for SvK models
476 // m_Gcov_ori_L = m_Gcov_ori;
477 // m_Gcov_ori.block(0,0,2,2) -= z*z * m_ncov_ori.transpose()*m_ncov_ori;
478 }
479
480 //--------------------------------------------------------------------------------------------------------------------------------------
481
482 template <short_t dim, class T>
_evalStretch(const gsMatrix<T> & C) const483 std::pair<gsVector<T>,gsMatrix<T>> gsMaterialMatrixBaseDim<dim,T>::_evalStretch(const gsMatrix<T> & C) const
484 {
485 gsVector<T> stretches;
486 gsMatrix<T> stretchvec;
487 std::pair<gsVector<T>,gsMatrix<T>> result;
488 stretches.resize(3,1); stretches.setZero();
489 stretchvec.resize(3,3); stretchvec.setZero();
490
491 Eigen::SelfAdjointEigenSolver< gsMatrix<real_t>::Base > eigSolver;
492
493
494 GISMO_ENSURE(m_gcon_ori.cols()!=0,"Is the basis initialized?");
495
496 gsMatrix<T,3,3> B;
497 B.setZero();
498 for (index_t k = 0; k != 2; k++)
499 for (index_t l = 0; l != 2; l++)
500 B += C(k,l) * m_gcon_ori.col(k) * m_gcon_ori.col(l).transpose();
501
502 eigSolver.compute(B);
503
504 stretchvec.leftCols(2) = eigSolver.eigenvectors().rightCols(2);
505 stretchvec.col(2) = m_gcon_ori.col(2); // replace with: stretchvec.col(0).template head<3>().cross(stretchvec.col(1).template head<3>())
506 stretches.block(0,0,2,1) = eigSolver.eigenvalues().block(1,0,2,1); // the eigenvalues are a 3x1 matrix, so we need to use matrix block-operations
507
508 // m_stretches.at(2) = 1/m_J0_sq;
509 stretches.at(2) = C(2,2);
510
511 for (index_t k=0; k!=3; k++)
512 stretches.at(k) = math::sqrt(stretches.at(k));
513
514 result.first = stretches;
515 result.second = stretchvec;
516
517 // // DEBUGGING ONLY!
518 // gsMatrix<T> ones(3,1);
519 // ones.setOnes();
520 // gsDebugVar(m_stretchvec);
521 // gsDebugVar(result.first);
522
523 return result;
524 }
525
526 //--------------------------------------------------------------------------------------------------------------------------------------
527
528 template <short_t dim, class T>
_computeStretch(const gsMatrix<T> & C) const529 void gsMaterialMatrixBaseDim<dim,T>::_computeStretch(const gsMatrix<T> & C) const
530 {
531 std::pair<gsVector<T>,gsMatrix<T>> result = _evalStretch(C);
532 m_stretches = result.first;
533 m_stretchvec = result.second;
534 }
535
536 //--------------------------------------------------------------------------------------------------------------------------------------
537
538 template <short_t dim, class T>
_transformation(const gsMatrix<T> & basis1,const gsMatrix<T> & basis2) const539 gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_transformation(const gsMatrix<T> & basis1, const gsMatrix<T> & basis2) const
540 {
541 // Transformation of a quantity FROM basis2 TO basis1
542 gsMatrix<T> Tmat(3,3);
543
544 for (index_t i = 0; i!=2; i++)
545 for (index_t j = 0; j!=2; j++)
546 Tmat(i,j) = math::pow(basis2.col(i).dot(basis1.col(j)),2);
547
548 Tmat(2,0) = basis2.col(1).dot(basis1.col(0)) * basis2.col(0).dot(basis1.col(0));
549 Tmat(2,1) = basis2.col(0).dot(basis1.col(1)) * basis2.col(1).dot(basis1.col(1));
550
551 Tmat(0,2) = 2*basis2.col(0).dot(basis1.col(1)) * basis2.col(0).dot(basis1.col(0));
552 Tmat(1,2) = 2*basis2.col(1).dot(basis1.col(0)) * basis2.col(1).dot(basis1.col(1));
553
554 Tmat(2,2) = basis2.col(0).dot(basis1.col(1)) * basis2.col(1).dot(basis1.col(0))
555 +basis2.col(0).dot(basis1.col(0)) * basis2.col(1).dot(basis1.col(1));
556
557 return Tmat;
558 }
559
560 } // end namespace
561