G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsMaterialMatrixComposite.hpp
Go to the documentation of this file.
1 
16 /*
17  To Do [updated 16-06-2020]:(rh)
18  - Make beta (compressible materials) and material parameters universal for all integration points over the thickness. So get them out of the dPsi functions etc and move them into the integration loops as global variables.
19 
20 */
21 
22 
23 
24 #pragma once
25 
26 #include <gsKLShell/src/gsMaterialMatrixComposite.h>
27 #include <gsKLShell/src/gsMaterialMatrixUtils.h>
28 #include <gsCore/gsFunction.h>
29 
30 namespace gismo
31 {
32 
33 template <short_t dim, class T >
34 gsMaterialMatrixComposite<dim,T>::gsMaterialMatrixComposite(
35  const gsFunctionSet<T> & mp,
36  const std::vector< gsFunctionSet<T> *> & thickness,
37  const std::vector< gsFunctionSet<T> *> & G,
38  const std::vector< gsFunctionSet<T> *> & alpha,
39  const std::vector< gsFunctionSet<T> *> & rho )
40  :
41  Base(&mp)
42 {
43  index_t nLayers = thickness.size();
44  GISMO_ASSERT(nLayers==(index_t)G.size(), "Size error in layer input");
45  GISMO_ASSERT(nLayers==(index_t)alpha.size(),"Size error in layer input");
46  GISMO_ASSERT(nLayers==(index_t)rho.size(), "Size error in layer input");
47 
48  for (size_t k=0; k!=thickness.size(); k++)
49  {
50  GISMO_ENSURE(thickness[k]!=NULL,"Function "<<k<<" of the thickness vector is not defined!");
51  GISMO_ENSURE(G[k]!=NULL,"Function "<<k<<" of the shear moduli vector is not defined!");
52  GISMO_ENSURE(alpha[k]!=NULL,"Function "<<k<<" of the alpha vector is not defined!");
53  }
54  for (size_t k=0; k!=rho.size(); k++)
55  GISMO_ENSURE(rho[k]!=NULL,"Function "<<k<<" of the density vector is not defined!");
56 
57  _initialize(nLayers);
58 
59  for (index_t k=0; k!=nLayers; k++)
60  {
61  m_Ts[k] = memory::make_shared_not_owned(thickness[k]);
62  m_Gs[k] = memory::make_shared_not_owned(G[k]);
63  m_As[k] = memory::make_shared_not_owned(alpha[k]);
64  m_Rs[k] = memory::make_shared_not_owned(rho[k]);
65  }
66 }
67 
68 template <short_t dim, class T >
69 gsMaterialMatrixComposite<dim,T>::gsMaterialMatrixComposite(
70  const gsFunctionSet<T> & mp,
71  const std::vector< gsFunctionSet<T> *> & thickness,
72  const std::vector< gsFunctionSet<T> *> & G,
73  const std::vector< gsFunctionSet<T> *> & alpha )
74  :
75  Base(&mp)
76 {
77  index_t nLayers = thickness.size();
78  GISMO_ASSERT(nLayers==(index_t)G.size(), "Size error in layer input");
79  GISMO_ASSERT(nLayers==(index_t)alpha.size(),"Size error in layer input");
80 
81  for (size_t k=0; k!=thickness.size(); k++)
82  {
83  GISMO_ENSURE(thickness[k]!=NULL,"Function "<<k<<" of the thickness vector is not defined!");
84  GISMO_ENSURE(G[k]!=NULL,"Function "<<k<<" of the shear moduli vector is not defined!");
85  GISMO_ENSURE(alpha[k]!=NULL,"Function "<<k<<" of the alpha vector is not defined!");
86  }
87 
88  _initialize(nLayers);
89  for (index_t k=0; k!=nLayers; k++)
90  {
91  m_Ts[k] = memory::make_shared_not_owned(thickness[k]);
92  m_Gs[k] = memory::make_shared_not_owned(G[k]);
93  m_As[k] = memory::make_shared_not_owned(alpha[k]);
94  }
95 }
96 
97 template <short_t dim, class T >
98 void gsMaterialMatrixComposite<dim,T>::_defaultOptions()
99 {
100  m_options.addInt("NumGauss","Number of Gaussian points through thickness",4);
101 }
102 
103 template <short_t dim, class T >
104 std::ostream & gsMaterialMatrixComposite<dim,T>::print(std::ostream &os) const
105 {
106  os <<"---------------------------------------------------------------------\n"
107  <<"----------------------Composite Material Info------------------------\n"
108  <<"---------------------------------------------------------------------\n\n";
109 
110  os <<"Material model: \t";
111  os <<"Saint-Venant Kirchhoff";
112  os <<"\n";
113 
114  os <<"---------------------------------------------------------------------\n\n";
115  return os;
116 }
117 
118 template <short_t dim, class T >
119 void gsMaterialMatrixComposite<dim,T>::_initialize(const index_t nLayers)
120 {
121  // Set default options
122  this->_defaultOptions();
123 
124  m_nLayers = nLayers;
125  m_Ts.resize(m_nLayers);
126  m_Gs.resize(m_nLayers);
127  m_As.resize(m_nLayers);
128  m_Rs.resize(m_nLayers);
129 }
130 
131 template <short_t dim, class T >
132 void gsMaterialMatrixComposite<dim,T>::_computePoints(const index_t patch, const gsMatrix<T> & u, bool basis) const
133 {
134  this->_computeMetricUndeformed(patch,u);
135 
136  if (Base::m_defpatches->nPieces()!=0)
137  this->_computeMetricDeformed(patch,u);
138 }
139 
140 template <short_t dim, class T >
141 void gsMaterialMatrixComposite<dim,T>::density_into(const index_t patch, const gsMatrix<T>& u, gsMatrix<T>& result) const
142 {
143  GISMO_ASSERT(m_Ts.size()==m_Rs.size(),"Size of vectors of thickness and densities is not equal: " << m_Ts.size()<<" & "<<m_Rs.size());
144 
145  gsMapData<T> map;
146  map.points = u;
147  map.flags = NEED_VALUE;
148  static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(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()
149 
150  result.resize(1, u.cols());
151  result.setZero();
152  std::vector<gsMatrix<T>> Tcontainer(m_nLayers);
153  std::vector<gsMatrix<T>> Rcontainer(m_nLayers);
154  for (size_t k=0; k!= m_Rs.size(); k++)
155  {
156  Tcontainer[k] = m_Ts[k]->piece(patch).eval(map.values[0]);
157  Rcontainer[k] = m_Rs[k]->piece(patch).eval(map.values[0]);
158 
159  GISMO_ASSERT(Tcontainer[k].rows()==1,"Thickness has the wrong size, must be scalar");
160  GISMO_ASSERT(Rcontainer[k].rows()==1,"Densities has the wrong size, must be scalar");
161  }
162 
163 
164  for (index_t i = 0; i != m_nLayers; ++i) // layers
165  for (index_t k = 0; k != u.cols(); ++k) // points
166  result(0,k) += Tcontainer[i](0,k)*Rcontainer[i](0,k);
167 
168 }
169 
170 template <short_t dim, class T >
171 void gsMaterialMatrixComposite<dim,T>::thickness_into(const index_t patch, const gsMatrix<T> & u, gsMatrix<T> & result) const
172 {
173  gsMapData<T> map;
174  map.points = u;
175  map.flags = NEED_VALUE;
176  static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(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()
177 
178  std::vector<gsMatrix<T>> Tcontainer(m_nLayers);
179  for (size_t k=0; k!= m_Rs.size(); k++)
180  {
181  Tcontainer[k] = m_Ts[k]->eval(map.values[0]);
182  GISMO_ASSERT(Tcontainer[k].rows()==1,"Thickness has the wrong size, must be scalar");
183  }
184 
185  for (index_t i = 0; i != m_nLayers; ++i) // layers
186  for (index_t k = 0; k != u.cols(); ++k) // points
187  result.row(i) = Tcontainer[i].row(0);
188 }
189 
190 template <short_t dim, class T >
191 gsMatrix<T> gsMaterialMatrixComposite<dim,T>::eval3D_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
192 {
193  // Input: u in-plane points
194  // z matrix with, per point, a column with z integration points
195  // Output: (n=u.cols(), m=z.cols())
196  // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
197 
198  this->_computePoints(patch,u);
199  // Initialize and result
200  gsMatrix<T> result(9, u.cols());
201  gsMatrix<T,3,3> Dmat, Cmat, Tmat;
202 
203  T t, t_tot, z_mid, t_temp, dz;
204 
205  std::vector<gsMatrix<T>> Tcontainer(m_nLayers);
206  std::vector<gsMatrix<T>> Rcontainer(m_nLayers);
207  std::vector<gsMatrix<T>> Acontainer(m_nLayers);
208  std::vector<gsMatrix<T>> Gcontainer(m_nLayers);
209 
210  gsMapData<T> map;
211  map.points = u;
212  map.flags = NEED_VALUE;
213  static_cast<const gsFunction<T>&>(m_patches->piece(patch) ).computeMap(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()
214 
215  gsMatrix<T> angles;
216  for (size_t k=0; k!= m_Gs.size(); k++)
217  {
218  Gcontainer[k] = m_Gs[k]->piece(patch).eval(map.values[0]);
219  Tcontainer[k] = m_Ts[k]->piece(patch).eval(map.values[0]);
220  angles = m_As[k]->piece(patch).eval(map.values[0]);
221 
222  GISMO_ASSERT(Gcontainer[k].rows()==9,"G has the wrong size, must be 3x3");
223  GISMO_ASSERT(Tcontainer[k].rows()==1,"Thickness has the wrong size, must be scalar");
224 
225  Acontainer[k] = _transformationMatrix(angles,map.values[0]);
226  }
227 
228 
229  for (index_t k = 0; k != u.cols(); ++k)
230  {
231  Cmat.setZero();
232  this->_getMetric(k,0.0); // on point i, with height 0.0
233 
234  // Compute total thickness (sum of entries)
235  t_tot = 0;
236  for (size_t i=0; i!=Tcontainer.size(); i++)
237  t_tot += Tcontainer[i](0,k);
238 
239  // compute mid-plane height of total plate
240  z_mid = t_tot / 2.0;
241 
242  // now we use t_temp to add the thickness of all plies iteratively
243  t_temp = 0.0;
244  t = 0.0;
245 
246  for (index_t i = 0; i != m_nLayers; ++i) // loop over laminates
247  {
248  // Lookup all quantities
249  t = Tcontainer[i](0,k);
250 
251  // Transform the matrix: T^T D T
252  Dmat = Acontainer[i].reshapeCol(k,3,3).transpose() * Gcontainer[i].reshapeCol(k,3,3) * Acontainer[i].reshapeCol(k,3,3);
253 
254  // distance from mid of the ply to mid of the plate
255  // dz = math::abs(z_mid - (t/2.0 + t_temp) ); // distance from mid-plane of plate
256  dz = z_mid - (t/2.0 + t_temp); // distance from mid-plane of plate
257 
258  if (out==MaterialOutput::MatrixA)
259  Cmat += Dmat * t; // A
260  else if (out==MaterialOutput::MatrixB || out==MaterialOutput::MatrixC)
261  Cmat += Dmat * t*dz; // B
262  else if (out==MaterialOutput::MatrixD)
263  Cmat += Dmat * ( t*math::pow(dz,2) + math::pow(t,3)/12.0 ); // D
264  else
265  GISMO_ERROR("MaterialOutput unknown");
266  t_temp += t;
267  }
268 
269  GISMO_ASSERT(std::abs(t_tot-t_temp)/t_tot < 1e-12,"Total thickness after loop is wrong. t_temp = "<<t_temp<<", sum(thickness) = "<<t_tot<<", difference = "<<t_tot-t_temp);
270 
271  gsVector<T,3> a1,ac1,e1,a2,ac2,e2;
272  a1 = m_data.mine().m_gcov_ori.col(0); a2 = m_data.mine().m_gcov_ori.col(1);
273  ac1 = m_data.mine().m_gcon_ori.col(0); ac2 = m_data.mine().m_gcon_ori.col(1);
274  e1 = m_data.mine().m_gcov_ori.col(0).normalized(); e2 = m_data.mine().m_gcon_ori.col(1).normalized();
275 
276  Cmat = _cart2cov(a1,a2,e1,e2) * Cmat;
277  Cmat = Cmat * _con2cart(ac1,ac2,e1,e2);
278 
279  result.reshapeCol(k,3,3) = Cmat;
280 
281  }
282  return result;
283 }
284 
285 template <short_t dim, class T >
286 gsMatrix<T> gsMaterialMatrixComposite<dim,T>::eval3D_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
287 {
288  // gsInfo<<"TO DO: evaluate moments using thickness";
289  // Input: u in-plane points
290  // z matrix with, per point, a column with z integration points
291  // Output: (n=u.cols(), m=z.cols())
292  // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
293 
294  this->_computePoints(patch,u);
295  gsMatrix<T> result(3, u.cols());
296  enum MaterialOutput _out;
297  if (out==MaterialOutput::VectorN)
298  _out=MaterialOutput::MatrixA;
299  else if (out==MaterialOutput::VectorM)
300  _out=MaterialOutput::MatrixD;
301  else
302  GISMO_ERROR("Output unknown");
303 
304  gsMatrix<T> Cmat = this->eval3D_matrix(patch,u,z,_out);
305  gsMatrix<T> Eij(2,2);
306 
307  for (index_t k = 0; k != u.cols(); ++k)
308  {
309  this->_getMetric(k,0.0); // on point i, with height 0.0
310 
311  if (out == MaterialOutput::VectorN) // To be used with multiplyZ_into
312  Eij = 0.5*(m_data.mine().m_Acov_def - m_data.mine().m_Acov_ori);
313  else if (out == MaterialOutput::VectorM) // To be used with multiplyZ_into
314  Eij = (m_data.mine().m_Bcov_ori - m_data.mine().m_Bcov_def);
315  else if (out == MaterialOutput::Generic) // To be used with multiplyLinZ_into or integrateZ_into
316  Eij = 0.5*(m_data.mine().m_Acov_def - m_data.mine().m_Acov_ori) + z*(m_data.mine().m_Bcov_ori - m_data.mine().m_Bcov_def);
317  else
318  GISMO_ERROR("Output type is not VectorN, VectorM or Generic!");
319 
320  Eij(0,1) += Eij(1,0);
321  std::swap(Eij(1,0), Eij(1,1));
322  Eij.resize(4,1);
323  Eij.conservativeResize(3,1);
324  result.col(k) = Cmat.reshapeCol(k,3,3) * Eij;
325  }
326  return result;
327 }
328 
329 //-----------------------------------------------------------------------------------------------------------------------
330 template <short_t dim, class T>
331 gsMatrix<T> gsMaterialMatrixComposite<dim,T>::_transformationMatrix(const gsMatrix<T> & phi, const gsMatrix<T> & u) const
332 {
333  gsMatrix<T> result(9,u.cols());
334  GISMO_ASSERT(phi.rows()==1,"Angles has the wrong size, must be scalar");
335 
336  for (index_t k=0; k!=u.cols(); k++)
337  {
338  gsAsMatrix<T,Dynamic,Dynamic> Tmat = result.reshapeCol(k,3,3);
339  // Make transformation matrix
340  Tmat(0,0) = Tmat(1,1) = math::pow(math::cos(phi(0,k)),2);
341  Tmat(0,1) = Tmat(1,0) = math::pow(math::sin(phi(0,k)),2);
342  Tmat(2,0) = Tmat(0,2) = Tmat(2,1) = Tmat(1,2) = math::sin(phi(0,k)) * math::cos(phi(0,k));
343  Tmat(2,0) *= -2.0;
344  Tmat(2,1) *= 2.0;
345  Tmat(1,2) *= -1.0;
346  Tmat(2,2) = math::pow(math::cos(phi(0,k)),2) - math::pow(math::sin(phi(0,k)),2);
347  }
348 
349  // Compute laminate stiffness matrix
350  return result;
351 }
352 
353 template <short_t dim, class T>
354 gsMatrix<T> gsMaterialMatrixComposite<dim,T>::_cart2cov(const gsVector<T> & a1, const gsVector<T> & a2, const gsVector<T> & e1, const gsVector<T> & e2) const
355 {
356  gsMatrix<T,3,3> Tmat;
357  Tmat.setZero();
358  // Covariant to local cartesian
359  Tmat(0,0) = (e1.dot(a1))*(a1.dot(e1));
360  Tmat(0,1) = (e1.dot(a2))*(a2.dot(e2));
361  Tmat(0,2) = 2*(e1.dot(a1))*(a2.dot(e1));
362  // Row 1
363  Tmat(1,0) = (e2.dot(a1))*(a1.dot(e2));
364  Tmat(1,1) = (e2.dot(a2))*(a2.dot(e2));
365  Tmat(1,2) = 2*(e2.dot(a1))*(a2.dot(e2));
366  // Row 2
367  Tmat(2,0) = (e1.dot(a1))*(a1.dot(e2));
368  Tmat(2,1) = (e1.dot(a2))*(a2.dot(e2));
369  Tmat(2,2) = (e1.dot(a1))*(a2.dot(e2)) + (e1.dot(a2))*(a1.dot(e2));
370 
371  Tmat = Tmat.template block<3,3>(0,0).inverse(); // !!!!
372  return Tmat;
373 }
374 
375 template <short_t dim, class T>
376 gsMatrix<T> gsMaterialMatrixComposite<dim,T>::_con2cart(const gsVector<T> & ac1, const gsVector<T> & ac2, const gsVector<T> & e1, const gsVector<T> & e2) const
377 {
378  gsMatrix<T,3,3> Tmat;
379  Tmat.setZero();
380  // Contravariant to local cartesian
381  Tmat(0,0) = (e1.dot(ac1))*(ac1.dot(e1));
382  Tmat(0,1) = (e1.dot(ac2))*(ac2.dot(e2));
383  Tmat(0,2) = 2*(e1.dot(ac1))*(ac2.dot(e1));
384  // Row 1
385  Tmat(1,0) = (e2.dot(ac1))*(ac1.dot(e2));
386  Tmat(1,1) = (e2.dot(ac2))*(ac2.dot(e2));
387  Tmat(1,2) = 2*(e2.dot(ac1))*(ac2.dot(e2));
388  // Row 2
389  Tmat(2,0) = (e1.dot(ac1))*(ac1.dot(e2));
390  Tmat(2,1) = (e1.dot(ac2))*(ac2.dot(e2));
391  Tmat(2,2) = (e1.dot(ac1))*(ac2.dot(e2)) + (e1.dot(ac2))*(ac1.dot(e2));
392  return Tmat;
393 }
394 
395 
396 } // end namespace
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition: gsMaterialMatrixUtils.h:98
gismo::memory::make_shared_not_owned
shared_ptr< T > make_shared_not_owned(const T *x)
Creates a shared pointer which does not eventually delete the underlying raw pointer. Usefull to refer to objects which should not be destroyed.
Definition: gsMemory.h:189
gismo::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition: gsLinearAlgebra.h:126
gismo::gsMapData
the gsMapData is a cache of pre-computed function (map) values.
Definition: gsFuncData.h:324
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition: gsFunction.h:59
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition: gsVector.h:35
gismo::gsMaterialMatrixComposite::thickness_into
void thickness_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixComposite.hpp:171
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
gsMaterialMatrixUtils.h
Provides material matrix utilities.
gismo::gsMaterialMatrixComposite::gsMaterialMatrixComposite
gsMaterialMatrixComposite()
Constructor of the assembler object.
Definition: gsMaterialMatrixComposite.h:56
gismo::gsMaterialMatrixComposite::eval3D_matrix
gsMatrix< T > eval3D_matrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
Evaluates the matrix on patch on in-plane points u with height z.
Definition: gsMaterialMatrixComposite.hpp:191
gismo::gsMaterialMatrixComposite::density_into
void density_into(const index_t patch, const gsMatrix< T > &u, gsMatrix< T > &result) const override
See gsMaterialMatrixBase for details.
Definition: gsMaterialMatrixComposite.hpp:141
gsMaterialMatrixComposite.h
Provides a material matrix for laminates.
gismo::NEED_VALUE
Value of the object.
Definition: gsForwardDeclarations.h:72
gismo::gsMaterialMatrixComposite
This class defines a linear material laminate.
Definition: gsMaterialMatrixComposite.h:39
GISMO_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
gismo::gsMapData::points
gsMatrix< T > points
input (parametric) points
Definition: gsFuncData.h:348
gismo::gsMaterialMatrixComposite::print
std::ostream & print(std::ostream &os) const override
Prints the object as a string.
Definition: gsMaterialMatrixComposite.hpp:104
gsFunction.h
Provides declaration of Function abstract interface.
gismo::gsMaterialMatrixComposite::eval3D_vector
gsMatrix< T > eval3D_vector(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
Evaluates the vector on patch on in-plane points u with height z.
Definition: gsMaterialMatrixComposite.hpp:286