G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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 <gsCore/gsGeometry.h>
27
28namespace gismo
29{
30
31//--------------------------------------------------------------------------------------------------------------------------------------
32
33template <short_t dim, class T >
34void gsMaterialMatrixBaseDim<dim,T>::defaultOptions()
35{
36 m_options.addInt("TensionField","Tension field definition - 0: Strain-based, 1: Stress-based, 2: Mixed",2);
37}
38
39template <short_t dim, class T >
40void gsMaterialMatrixBaseDim<dim,T>::density_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
41{
42 gsMapData<T> map;
43 map.flags = NEED_VALUE;
44 map.points = u;
45 static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
46
47 result.resize(1, u.cols());
48 m_thickness->piece(patch).eval_into(map.values[0], m_data.mine().m_Tmat);
49 m_density->piece(patch).eval_into(map.values[0], m_data.mine().m_rhomat);
50 for (index_t i = 0; i != u.cols(); ++i) // points
51 {
52 result(0,i) = m_data.mine().m_Tmat(0,i)*m_data.mine().m_rhomat(0,i);
53 }
54}
55
56template <short_t dim, class T >
57void gsMaterialMatrixBaseDim<dim,T>::thickness_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
58{
59 gsMapData<T> map;
60 map.flags = NEED_VALUE;
61 map.points = u;
62 static_cast<const gsFunction<T>&>(Base::m_patches->piece(patch) ).computeMap(map);
63 m_thickness->piece(patch).eval_into(map.values[0], result);
64}
65
66// Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) covariant basis
67template <short_t dim, class T>
68gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_spec2cov(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
69{
70 gsMatrix<T> result(9, u.cols());
71 gsMatrix<T> covbasis,sbasis;
72 gsMatrix<T> tmp = eval3D_pstretchDir(patch,u,z);
73 index_t colIdx;
74 for (index_t k=0; k!=u.cols(); k++)
75 {
76 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
77 {
78 colIdx = j*u.cols()+k;
79 _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
80 sbasis = tmp.reshapeCol(colIdx,3,3);
81 covbasis = m_data.mine().m_gcov_ori;
82 result.col(colIdx) = this->_transformation(covbasis,sbasis).reshape(9,1);
83 }
84 }
85 return result;
86}
87
88// Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) convariant basis
89template <short_t dim, class T>
90gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_spec2con(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
91{
92 gsMatrix<T> result(9, u.cols());
93 gsMatrix<T> conbasis,sbasis;
94 gsMatrix<T> tmp = eval3D_pstretchDir(patch,u,z);
95 index_t colIdx;
96 for (index_t k=0; k!=u.cols(); k++)
97 {
98 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
99 {
100 colIdx = j*u.cols()+k;
101 _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
102 sbasis = tmp.reshapeCol(colIdx,3,3);
103 conbasis = m_data.mine().m_gcov_ori;
104 result.col(colIdx) = this->_transformation(conbasis,sbasis).reshape(9,1);
105 }
106 }
107 return result;
108}
109
110// Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) covariant basis
111template <short_t dim, class T>
112gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_cov2cart(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
113{
114 this->_computePoints(patch,u);
115
116 gsMatrix<T> result(9, u.cols());
117 gsMatrix<T> tmp, covbasis,cartbasis(3,3);
118 cartbasis.setIdentity();
119 index_t colIdx;
120 for (index_t k=0; k!=u.cols(); k++)
121 {
122 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
123 {
124 colIdx = j*u.cols()+k;
125 _getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
126 covbasis = m_data.mine().m_gcov_ori;
127 result.col(colIdx) = this->_transformation(cartbasis,covbasis).reshape(9,1);
128 }
129 }
130 return result;
131}
132
133// Constructs a transformation matrix that transforms a quantity (IN VOIGHT NOTATION) in the spectral basis to the (undeformed) convariant basis
134template <short_t dim, class T>
135gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::eval3D_con2cart(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
136{
137 this->_computePoints(patch,u);
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
156template <short_t dim, class T >
157void 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
175template <short_t dim, class T >
176gsMatrix<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
199template <short_t dim, class T>
200gsMatrix<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
221template <short_t dim, class T>
222gsMatrix<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
243template <short_t dim, class T>
244gsMatrix<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
270template <short_t dim, class T>
271gsMatrix<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
295template <short_t dim, class T>
296gsMatrix<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
345template <short_t dim, class T>
346template <short_t _type>
347typename std::enable_if<_type==0, index_t>::type
348gsMaterialMatrixBaseDim<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
358template <short_t dim, class T>
359template <short_t _type>
360typename std::enable_if<_type==1, index_t>::type
361gsMaterialMatrixBaseDim<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
371template <short_t dim, class T>
372template <short_t _type>
373typename std::enable_if<_type==2, index_t>::type
374gsMaterialMatrixBaseDim<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
388template <short_t dim, class T >
389void 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
416template <short_t dim, class T>
417void gsMaterialMatrixBaseDim<dim,T>::_computeMetricDeformed(const index_t patch, const gsMatrix<T> & u) const
418{
419 _computeMetricDeformed_impl<dim>(patch,u);
420}
421
422template <short_t dim, class T>
423template <short_t _dim>
424typename std::enable_if<_dim==3, void>::type
425gsMaterialMatrixBaseDim<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
485template <short_t dim, class T>
486template <short_t _dim>
487typename std::enable_if<_dim==2, void>::type
488gsMaterialMatrixBaseDim<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
534template <short_t dim, class T>
535void gsMaterialMatrixBaseDim<dim,T>::_computeMetricUndeformed(const index_t patch, const gsMatrix<T> & u) const
536{
537 _computeMetricUndeformed_impl<dim>(patch,u);
538}
539
540template <short_t dim, class T>
541template <short_t _dim>
542typename std::enable_if<_dim==3, void>::type
543gsMaterialMatrixBaseDim<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
603template <short_t dim, class T>
604template <short_t _dim>
605typename std::enable_if<_dim==2, void>::type
606gsMaterialMatrixBaseDim<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
650template <short_t dim, class T>
651gsMatrix<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
657template <short_t dim, class T>
658gsMatrix<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
664template <short_t dim, class T>
665gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getBcov_def(index_t k, T z) const
666{
667 return _getBcov_def_impl<dim>(k,z);
668}
669
670template <short_t dim, class T>
671gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getncov_def(index_t k, T z) const
672{
673 return _getncov_def_impl<dim>(k,z);
674}
675
676template <short_t dim, class T>
677gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcov_def(index_t k, T z) const
678{
679 return _getGcov_def_impl<dim>(k,z);
680}
681
682template <short_t dim, class T>
683gsMatrix<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
689template <short_t dim, class T>
690gsMatrix<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
696template <short_t dim, class T>
697gsMatrix<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
703template <short_t dim, class T>
704gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcov_def(index_t k, T z) const
705{
706 return _getgcov_def_impl<dim>(k,z);
707}
708
709template <short_t dim, class T>
710gsMatrix<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
724template <short_t dim, class T>
725gsMatrix<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
731template <short_t dim, class T>
732gsMatrix<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
738template <short_t dim, class T>
739gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getBcov_ori(index_t k, T z) const
740{
741 return _getBcov_ori_impl<dim>(k,z);
742}
743
744template <short_t dim, class T>
745gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getncov_ori(index_t k, T z) const
746{
747 return _getncov_ori_impl<dim>(k,z);
748}
749
750template <short_t dim, class T>
751gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getGcov_ori(index_t k, T z) const
752{
753 return _getGcov_ori_impl<dim>(k,z);
754}
755
756template <short_t dim, class T>
757gsMatrix<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
763template <short_t dim, class T>
764gsMatrix<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
770template <short_t dim, class T>
771gsMatrix<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
777template <short_t dim, class T>
778gsMatrix<T> gsMaterialMatrixBaseDim<dim,T>::_getgcov_ori(index_t k, T z) const
779{
780 return _getgcov_ori_impl<dim>(k,z);
781}
782
783template <short_t dim, class T>
784gsMatrix<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
800template <short_t dim, class T>
801template <short_t _dim>
802typename std::enable_if<_dim==2, gsMatrix<T>>::type
807
808template <short_t dim, class T>
809template <short_t _dim>
810typename std::enable_if<_dim==3, gsMatrix<T>>::type
811gsMaterialMatrixBaseDim<dim,T>::_getBcov_def_impl(index_t k, T) 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
817template <short_t dim, class T>
818template <short_t _dim>
819typename std::enable_if<_dim==2, gsMatrix<T>>::type
824
825template <short_t dim, class T>
826template <short_t _dim>
827typename std::enable_if<_dim==3, gsMatrix<T>>::type
828gsMaterialMatrixBaseDim<dim,T>::_getncov_def_impl(index_t k, T) 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
834template <short_t dim, class T>
835template <short_t _dim>
836typename std::enable_if<_dim==2, gsMatrix<T>>::type
837gsMaterialMatrixBaseDim<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
850template <short_t dim, class T>
851template <short_t _dim>
852typename std::enable_if<_dim==3, gsMatrix<T>>::type
853gsMaterialMatrixBaseDim<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
868template <short_t dim, class T>
869template <short_t _dim>
870typename std::enable_if<_dim==2, gsMatrix<T>>::type
871gsMaterialMatrixBaseDim<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
882template <short_t dim, class T>
883template <short_t _dim>
884typename std::enable_if<_dim==3, gsMatrix<T>>::type
885gsMaterialMatrixBaseDim<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
898template <short_t dim, class T>
899template <short_t _dim>
900typename std::enable_if<_dim==2, gsMatrix<T>>::type
905
906template <short_t dim, class T>
907template <short_t _dim>
908typename std::enable_if<_dim==3, gsMatrix<T>>::type
909gsMaterialMatrixBaseDim<dim,T>::_getBcov_ori_impl(index_t k, T) 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
915template <short_t dim, class T>
916template <short_t _dim>
917typename std::enable_if<_dim==2, gsMatrix<T>>::type
922
923template <short_t dim, class T>
924template <short_t _dim>
925typename std::enable_if<_dim==3, gsMatrix<T>>::type
926gsMaterialMatrixBaseDim<dim,T>::_getncov_ori_impl(index_t k, T) 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
932template <short_t dim, class T>
933template <short_t _dim>
934typename std::enable_if<_dim==2, gsMatrix<T>>::type
935gsMaterialMatrixBaseDim<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
948template <short_t dim, class T>
949template <short_t _dim>
950typename std::enable_if<_dim==3, gsMatrix<T>>::type
951gsMaterialMatrixBaseDim<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
966template <short_t dim, class T>
967template <short_t _dim>
968typename std::enable_if<_dim==2, gsMatrix<T>>::type
969gsMaterialMatrixBaseDim<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
980template <short_t dim, class T>
981template <short_t _dim>
982typename std::enable_if<_dim==3, gsMatrix<T>>::type
983gsMaterialMatrixBaseDim<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
998template <short_t dim, class T>
999void 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
1020template <short_t dim, class T>
1021void 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
1044template <short_t dim, class T>
1045void 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
1072template <short_t dim, class T>
1073void gsMaterialMatrixBaseDim<dim,T>::_getMetricDeformed(const index_t k, const T z) const
1074{
1075 _getMetricDeformed_impl<dim>(k,z);
1076}
1077
1078template <short_t dim, class T>
1079template <short_t _dim>
1080typename std::enable_if<_dim==3, void>::type
1081gsMaterialMatrixBaseDim<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
1138template <short_t dim, class T>
1139template <short_t _dim>
1140typename std::enable_if<_dim==2, void>::type
1141gsMaterialMatrixBaseDim<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
1201template <short_t dim, class T>
1202void gsMaterialMatrixBaseDim<dim,T>::_getMetricUndeformed(const index_t k, const T z) const
1203{
1204 _getMetricUndeformed_impl<dim>(k,z);
1205}
1206
1207template <short_t dim, class T>
1208template <short_t _dim>
1209typename std::enable_if<_dim==3, void>::type
1210gsMaterialMatrixBaseDim<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
1269template <short_t dim, class T>
1270template <short_t _dim>
1271typename std::enable_if<_dim==2, void>::type
1272gsMaterialMatrixBaseDim<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
1330template <short_t dim, class T>
1331std::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
1383template <short_t dim, class T>
1384std::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
1427template <short_t dim, class T>
1428std::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
1475template <short_t dim, class T>
1476void 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
1483template <short_t dim, class T>
1484void 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
1491template <short_t dim, class T>
1492void 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
1501template <short_t dim, class T>
1502gsMatrix<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::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition gsAsMatrix.h:32
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition gsFunction.h:60
gismo::gsMapData
the gsMapData is a cache of pre-computed function (map) values.
Definition gsFuncData.h:349
gismo::gsMapData::points
gsMatrix< T > points
input (parametric) points
Definition gsFuncData.h:372
gismo::gsMaterialMatrixBaseDim
This class defines the base class for material matrices.
Definition gsMaterialMatrixBaseDim.h:37
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::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::_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::_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::_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::_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:543
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::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:90
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::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:425
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
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:40
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::_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::_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::_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:1081
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::_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_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:135
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::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::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::_computeMetricDeformed
void _computeMetricDeformed(const index_t patch, const gsMatrix< T > &u) const
Computes metric quantities on the deformed geometry.
Definition gsMaterialMatrixBaseDim.hpp:417
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::_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::_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::_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::defaultOptions
virtual void defaultOptions() override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixBaseDim.hpp:34
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::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::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::_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::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::_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
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 basis vector g on the deformed geometry.
Definition gsMaterialMatrixBaseDim.hpp:704
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::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::_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::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:57
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::_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::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:112
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::_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:1210
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::_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::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:68
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::_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::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::_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::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
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::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::_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::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshapeCol
gsAsMatrix< T, Dynamic, Dynamic > reshapeCol(index_t c, index_t n, index_t m)
Returns column c of the matrix resized to n x m matrix This function assumes that the matrix is size ...
Definition gsMatrix.h:231
gismo::gsMatrix::reshape
gsAsMatrix< T, Dynamic, Dynamic > reshape(index_t n, index_t m)
Returns the matrix resized to n x m matrix (data is not copied) This function assumes that the matrix...
Definition gsMatrix.h:221
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition gsVector.h:177
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition gsMaterialMatrixUtils.h:99
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsGeometry.h
Provides declaration of Geometry abstract interface.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::NEED_VALUE
@ NEED_VALUE
Value of the object.
Definition gsForwardDeclarations.h:72
gismo::NEED_DERIV2
@ NEED_DERIV2
Second derivatives.
Definition gsForwardDeclarations.h:80
gismo::NEED_DERIV
@ NEED_DERIV
Gradient of the object.
Definition gsForwardDeclarations.h:73
gismo::NEED_JACOBIAN
@ NEED_JACOBIAN
Jacobian of the object.
Definition gsForwardDeclarations.h:75
gismo::NEED_NORMAL
@ NEED_NORMAL
Normal vector of the object.
Definition gsForwardDeclarations.h:85