G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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
30namespace gismo
31{
32
33template <short_t dim, class T >
34gsMaterialMatrixComposite<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
68template <short_t dim, class T >
69gsMaterialMatrixComposite<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
97template <short_t dim, class T >
98void gsMaterialMatrixComposite<dim,T>::_defaultOptions()
99{
100 m_options.addInt("NumGauss","Number of Gaussian points through thickness",4);
101}
102
103template <short_t dim, class T >
104std::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
118template <short_t dim, class T >
119void 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
131template <short_t dim, class T >
132void 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
140template <short_t dim, class T >
141void 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
170template <short_t dim, class T >
171void 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
190template <short_t dim, class T >
191gsMatrix<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
285template <short_t dim, class T >
286gsMatrix<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//-----------------------------------------------------------------------------------------------------------------------
330template <short_t dim, class T>
331gsMatrix<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
353template <short_t dim, class T>
354gsMatrix<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
375template <short_t dim, class T>
376gsMatrix<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::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::gsMaterialMatrixComposite
This class defines a linear material laminate.
Definition gsMaterialMatrixComposite.h:40
gismo::gsMaterialMatrixComposite::gsMaterialMatrixComposite
gsMaterialMatrixComposite()
Constructor of the assembler object.
Definition gsMaterialMatrixComposite.h:56
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::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition gsMaterialMatrixUtils.h:99
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsFunction.h
Provides declaration of Function abstract interface.
gsMaterialMatrixComposite.h
Provides a material matrix for laminates.
gsMaterialMatrixUtils.h
Provides material matrix utilities.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::NEED_VALUE
@ NEED_VALUE
Value of the object.
Definition gsForwardDeclarations.h:72