G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMaterialMatrixIntegrate.hpp
Go to the documentation of this file.
1
16/*
17 To Do [updated 16-06-2020]:
18 - Make beta (compressible materials) and material parameters universal for all integration points over the thickness. So get them out of the dPsi functions etc and move them into the integration loops as global variables.
19
20*/
21
22
23
24#pragma once
25
26#include <gsKLShell/src/gsMaterialMatrixIntegrate.h>
27#include <gsAssembler/gsGaussRule.h>
28
29namespace gismo
30{
31
32template <class T, enum MaterialOutput out>
33gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(const gsMaterialMatrixContainer<T> & materialMatrices,
34 const gsFunctionSet<T> * deformed)
35:
36m_materialMatrices(materialMatrices),
37m_deformed(deformed)
38{
39 for (index_t p = 0; p!=m_materialMatrices.size(); p++)
40 {
41 GISMO_ASSERT(materialMatrices.piece(p)!=nullptr,"Material matrix "<<p<<" is incomplete!");
42 GISMO_ASSERT(materialMatrices.piece(p)!=NULL,"Material matrix "<<p<<" is incomplete!");
43 GISMO_ASSERT(materialMatrices.piece(p)->initialized(),"Material matrix "<<p<<" is incomplete!");
44 }
45
46 this->_makePieces(deformed);
47}
48
49template <class T, enum MaterialOutput out>
50gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(gsMaterialMatrixBase<T> * materialMatrix,
51 const gsFunctionSet<T> * deformed)
52:
53m_materialMatrices(deformed->nPieces()),
54m_deformed(deformed)
55{
56 GISMO_ASSERT(materialMatrix->initialized(),"Material matrix is incomplete!");
57 for (index_t p = 0; p!=deformed->nPieces(); ++p)
58 m_materialMatrices.set(p,materialMatrix);
59 this->_makePieces(deformed);
60}
61
62template <class T, enum MaterialOutput out>
63gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(typename gsMaterialMatrixBase<T>::uPtr & materialMatrix,
64 const gsFunctionSet<T> * deformed)
65:
66gsMaterialMatrixIntegrate(materialMatrix.get(),deformed)
67{
68}
69
70template <class T, enum MaterialOutput out>
71gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(const gsMaterialMatrixContainer<T> & materialMatrices,
72 const gsFunctionSet<T> * undeformed,
73 const gsFunctionSet<T> * deformed)
74:
75m_materialMatrices(materialMatrices),
76m_deformed(deformed)
77{
78 this->_makePieces(undeformed,deformed);
79}
80
81template <class T, enum MaterialOutput out>
82gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(gsMaterialMatrixBase<T> * materialMatrix,
83 const gsFunctionSet<T> * undeformed,
84 const gsFunctionSet<T> * deformed)
85:
86m_materialMatrices(deformed->nPieces()),
87m_deformed(deformed)
88{
89 for (index_t p = 0; p!=deformed->nPieces(); ++p)
90 {
91 m_materialMatrices.set(p,materialMatrix);
92 this->_makePieces(undeformed,deformed);
93 }
94}
95
96template <class T, enum MaterialOutput out>
97gsMaterialMatrixIntegrate<T,out>::gsMaterialMatrixIntegrate(typename gsMaterialMatrixBase<T>::uPtr & materialMatrix,
98 const gsFunctionSet<T> * undeformed,
99 const gsFunctionSet<T> * deformed)
100:
101gsMaterialMatrixIntegrate(materialMatrix.get(),undeformed,deformed)
102{
103}
104
105template <class T, enum MaterialOutput out>
106gsMaterialMatrixIntegrate<T,out>::~gsMaterialMatrixIntegrate()
107{
108 freeAll(m_pieces);
109}
110
111template <class T, enum MaterialOutput out>
112void gsMaterialMatrixIntegrate<T,out>::_makePieces(const gsFunctionSet<T> * deformed)
113{
114 m_pieces.resize(deformed->nPieces());
115 for (size_t p = 0; p!=m_pieces.size(); ++p)
116 m_pieces.at(p) = new gsMaterialMatrixIntegrateSingle<T,out>(p,m_materialMatrices.piece(p),deformed);
117}
118
119
120template <class T, enum MaterialOutput out>
121void gsMaterialMatrixIntegrate<T,out>::_makePieces(const gsFunctionSet<T> * undeformed, const gsFunctionSet<T> * deformed)
122{
123 m_pieces.resize(m_deformed->nPieces());
124 for (size_t p = 0; p!=m_pieces.size(); ++p)
125 m_pieces.at(p) = new gsMaterialMatrixIntegrateSingle<T,out>(p,m_materialMatrices.piece(p),undeformed,deformed);
126}
127
128
129// Linear material models
130template <class T, enum MaterialOutput out>
131gsMaterialMatrixIntegrateSingle<T,out>::gsMaterialMatrixIntegrateSingle(index_t patch,
132 gsMaterialMatrixBase<T> * materialMatrix, //??
133 const gsFunctionSet<T> * deformed
134 )
135:
136m_pIndex(patch),
137m_materialMat(materialMatrix)
138{
139 m_materialMat->setDeformed(deformed);
140 // m_materialMat = new gsMaterialMatrix(materialMatrix);
141}
142
143// Linear material models
144template <class T, enum MaterialOutput out>
145gsMaterialMatrixIntegrateSingle<T,out>::gsMaterialMatrixIntegrateSingle(index_t patch,
146 gsMaterialMatrixBase<T> * materialMatrix, //??
147 const gsFunctionSet<T> * undeformed,
148 const gsFunctionSet<T> * deformed
149 )
150:
151m_pIndex(patch),
152m_materialMat(materialMatrix)
153{
154 m_materialMat->setUndeformed(undeformed);
155 m_materialMat->setDeformed(deformed);
156 // m_materialMat = new gsMaterialMatrix(materialMatrix);
157}
158
159template <class T, enum MaterialOutput out>
160void gsMaterialMatrixIntegrateSingle<T,out>::eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
161{
163// #pragma omp critical (gsMaterialMatrixIntegrateSingle_eval_into)
164 this->eval_into_impl<out>(u,result);
165}
166
167template <class T, enum MaterialOutput out>
168template <enum MaterialOutput _out>
169typename std::enable_if<_out==MaterialOutput::Density, void>::type
170gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
171{
172 m_materialMat->density_into(m_pIndex,u,result);
173}
174
175template <class T, enum MaterialOutput out>
176template <enum MaterialOutput _out>
177typename std::enable_if<_out==MaterialOutput::VectorN ||
178 _out==MaterialOutput::CauchyVectorN, void>::type
179gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
180{
181 if (m_materialMat->isVecIntegrated() == MatIntegration::NotIntegrated)
182 this->integrateZ_into(u,getMoment(),result);
183 else if (m_materialMat->isVecIntegrated() == MatIntegration::Integrated)
184 result = this->_eval(u);
185 else if (m_materialMat->isVecIntegrated() == MatIntegration::Constant)
186 this->multiplyZ_into(u,0,result);
187 else if (m_materialMat->isVecIntegrated() == MatIntegration::Linear)
188 this->multiplyLinZ_into(u,getMoment(),result);
189 else
190 GISMO_ERROR("Integration status unknown");
191}
192
193template <class T, enum MaterialOutput out>
194template <enum MaterialOutput _out>
195typename std::enable_if<_out==MaterialOutput::VectorM ||
196 _out==MaterialOutput::CauchyVectorM, void>::type
197gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
198{
199 if (m_materialMat->isVecIntegrated() == MatIntegration::NotIntegrated)
200 this->integrateZ_into(u,getMoment(),result);
201 else if (m_materialMat->isVecIntegrated() == MatIntegration::Integrated)
202 result = this->_eval(u);
203 else if (m_materialMat->isVecIntegrated() == MatIntegration::Constant)
204 this->multiplyZ_into(u,2,result);
205 else if (m_materialMat->isVecIntegrated() == MatIntegration::Linear)
206 this->multiplyLinZ_into(u,getMoment(),result);
207 else
208 GISMO_ERROR("Integration status unknown");
209}
210
211template <class T, enum MaterialOutput out>
212template <enum MaterialOutput _out>
213typename std::enable_if< _out==MaterialOutput::MatrixA || _out==MaterialOutput::MatrixB
214 || _out==MaterialOutput::MatrixC || _out==MaterialOutput::MatrixD, void>::type
215gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
216{
217 if (m_materialMat->isMatIntegrated() == MatIntegration::NotIntegrated)
218 this->integrateZ_into(u,getMoment(),result);
219 else if (m_materialMat->isMatIntegrated() == MatIntegration::Integrated)
220 result = this->_eval(u);
221 else if (m_materialMat->isMatIntegrated() == MatIntegration::Constant)
222 this->multiplyZ_into(u,getMoment(),result);
223 else if (m_materialMat->isMatIntegrated() == MatIntegration::Linear)
224 this->multiplyLinZ_into(u,getMoment(),result);
225 else
226 GISMO_ERROR("Integration status unknown");
227}
228
229template <class T, enum MaterialOutput out>
230template <enum MaterialOutput _out>
231typename std::enable_if<_out==MaterialOutput::PStressN || _out==MaterialOutput::PStressM, void>::type
232gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
233{
234 if (m_materialMat->isMatIntegrated() == MatIntegration::NotIntegrated)
235 this->integrateZ_into(u,getMoment(),result);
236 else if (m_materialMat->isMatIntegrated() == MatIntegration::Integrated)
237 result = this->_eval(u);
238 else if (m_materialMat->isMatIntegrated() == MatIntegration::Constant)
239 this->multiplyZ_into(u,getMoment(),result);
240 else if (m_materialMat->isMatIntegrated() == MatIntegration::Linear)
241 this->multiplyLinZ_into(u,getMoment(),result);
242 else
243 GISMO_ERROR("Integration status unknown");
244}
245
246template <class T, enum MaterialOutput out>
247template <enum MaterialOutput _out>
248typename std::enable_if<_out==MaterialOutput::PStrainN || _out==MaterialOutput::PStrainM, void>::type
249gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
250{
251 if (m_materialMat->isMatIntegrated() == MatIntegration::NotIntegrated)
252 this->integrateZ_into(u,getMoment(),result);
253 else if (m_materialMat->isMatIntegrated() == MatIntegration::Integrated)
254 result = this->_eval(u);
255 else if (m_materialMat->isMatIntegrated() == MatIntegration::Constant)
256 this->multiplyZ_into(u,getMoment(),result);
257 else if (m_materialMat->isMatIntegrated() == MatIntegration::Linear)
258 this->multiplyLinZ_into(u,getMoment(),result);
259 else
260 GISMO_ERROR("Integration status unknown");
261}
262
263template <class T, enum MaterialOutput out>
264template <enum MaterialOutput _out>
265typename std::enable_if<_out==MaterialOutput::Stretch, void>::type
266gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
267{
268 m_materialMat->stretch_into(m_pIndex,u,result);
269}
270
271template <class T, enum MaterialOutput out>
272template <enum MaterialOutput _out>
273typename std::enable_if<_out==MaterialOutput::StretchDir, void>::type
274gsMaterialMatrixIntegrateSingle<T,out>::eval_into_impl(const gsMatrix<T>& u, gsMatrix<T>& result) const
275{
276 m_materialMat->stretchDir_into(m_pIndex,u,result);
277}
278
279template <class T, enum MaterialOutput out>
280void gsMaterialMatrixIntegrateSingle<T,out>::integrateZ_into(const gsMatrix<T>& u, const index_t moment, gsMatrix<T> & result) const
281{
282 // Input: points in R2
283 // Ouput: results in targetDim
284
285 // Perform integration
286 index_t numGauss = m_materialMat->options().askInt("NumGauss",4);
287
288 result.resize(this->targetDim(),u.cols());
289 result.setZero();
290
291 gsGaussRule<T> gauss(numGauss);
292 gsMatrix<T> z(numGauss,u.cols());
293 gsMatrix<T> w(numGauss,u.cols());
294 gsMatrix<T> vals;
295 // m_points3D.resize(1,numGauss);
296
297 gsMatrix<T> Tmat;
298 m_materialMat->thickness_into(m_pIndex,u,Tmat);
299 T Thalf;
300 // pre-compute Z
301 for (index_t k = 0; k != u.cols(); ++k) // for all points
302 {
303 gsMatrix<T> quNodes(1,numGauss);
304 gsVector<T> quWeights(numGauss);
305 // set new integration point
306 // Thalf = Tmat(0,k)/2.0;
307 Thalf = 0.5;
308 gauss.mapTo(-Thalf,Thalf,quNodes,quWeights);
309 w.col(k)=quWeights;
310 z.col(k)=quNodes.transpose();
311 }
312 vals = this->_eval3D(u,z);
313 T res;
314 for (index_t k = 0; k != u.cols(); ++k) // for all points
315 {
316 for (index_t i=0; i!=this->targetDim(); ++i) // components
317 {
318 res = 0.0;
319 for (index_t j = 0; j != numGauss; ++j) // compute integral
320 res += w(j, k) * math::pow(z(j, k) * Tmat(0, k), moment) * vals(i, j * u.cols() + k) * Tmat(0, k); // int_z=[-1/2,1/2] (xi*t)^moment * <quantity(xi)> * t
321 result(i,k) = res;
322 }
323
324 }
325}
326
327
328/*
329 * WARNING: This function assumes the function to be integrated: f(z,...) = g(,...) + z h(...)
330 Note: z is assumed to be [-1,1]
331 */
332template <class T, enum MaterialOutput out>
333void gsMaterialMatrixIntegrateSingle<T,out>::multiplyLinZ_into(const gsMatrix<T>& u, const index_t moment, gsMatrix<T>& result) const
334{
335 // Input: points in R2
336 // Ouput: results in targetDim
337 result.resize(this->targetDim(),u.cols());
338
339 gsMatrix<T> z(2,u.cols());
340 gsMatrix<T> Tmat;
341 m_materialMat->thickness_into(m_pIndex,u,Tmat);
342
343 T Thalf;
344 // pre-compute Z
345 for (index_t k = 0; k != u.cols(); ++k) // for all points
346 {
347 Thalf = 1.0 / 2.0;
348 // Thalf = Tmat(0, k) / 2.0;
349 z.col(k)<<Thalf,-Thalf;
350 }
351
352 T fac = (moment % 2 == 0) ? 0. : 1.;
353
354 gsMatrix<T> vals = this->_eval3D(u,z);
355 for (index_t k = 0; k != u.cols(); ++k) // for all points
356 {
357 // 1/(alpha+1) * [ (t/2)^(alpha+1) * g(...) - (-t/2)^(alpha+1) * g(...) ]
358 // 1/(alpha+2) * [ (t/2)^(alpha+1) * h(...) - (-t/2)^(alpha+1) * h(...) ]
359 result.col(k) = 1.0/(moment+fac+1) *
360 ( math::pow( z(0,k) * Tmat(0, k), moment + 1) * vals.col(0*u.cols() + k)
361 - math::pow( z(1,k) * Tmat(0, k), moment + 1) * vals.col(1*u.cols() + k) );
362 }
363
364}
365
366template <class T, enum MaterialOutput out>
367void gsMaterialMatrixIntegrateSingle<T,out>::multiplyZ_into(const gsMatrix<T>& u, index_t moment, gsMatrix<T>& result) const
368{
369 // Input: points in R2
370 // Ouput: results in targetDim
371 result.resize(this->targetDim(),u.cols());
372
373 if (moment % 2 != 0) //then the moment is odd
374 result.setZero();
375 else // then the moment is even
376 {
377 T fac;
378 gsMatrix<T> Tmat;
379 m_materialMat->thickness_into(m_pIndex,u,Tmat);
380 T Thalf;
381 gsMatrix<T> vals = this->_eval(u);
382 for (index_t k = 0; k != u.cols(); ++k) // for all points
383 {
384 Thalf = Tmat(0, k) / 2.0;
385 fac = 2.0/(moment+1) * math::pow( Thalf , moment + 1);
386 result.col(k) = vals.col(k) * fac;
387 }
388
389 }
390}
391
392template <class T, enum MaterialOutput out>
393gsMatrix<T> gsMaterialMatrixIntegrateSingle<T,out>::_eval(const gsMatrix<T>& u) const
394{
395 gsMatrix<T> Z(1,u.cols());
396 Z.setZero();
397 return this->_eval3D(u,Z);
398}
399
400template <class T, enum MaterialOutput out>
401gsMatrix<T> gsMaterialMatrixIntegrateSingle<T,out>::_eval3D(const gsMatrix<T>& u, const gsMatrix<T>& Z) const
402{
403 return this->eval3D_impl<out>(u,Z);
404}
405
406template <class T, enum MaterialOutput out>
407template <enum MaterialOutput _out>
408typename std::enable_if< _out==MaterialOutput::MatrixA || _out==MaterialOutput::MatrixB
409 || _out==MaterialOutput::MatrixC || _out==MaterialOutput::MatrixD, gsMatrix<T>>::type
410gsMaterialMatrixIntegrateSingle<T,out>::eval3D_impl(const gsMatrix<T>& u, const gsMatrix<T>& Z) const
411{
412 return m_materialMat->eval3D_matrix(m_pIndex,u,Z,_out);
413}
414
415template <class T, enum MaterialOutput out>
416template <enum MaterialOutput _out>
417typename std::enable_if<_out==MaterialOutput::VectorN ||
418 _out==MaterialOutput::VectorM , gsMatrix<T>>::type
419gsMaterialMatrixIntegrateSingle<T,out>::eval3D_impl(const gsMatrix<T>& u, const gsMatrix<T>& Z) const
420{
421 return m_materialMat->eval3D_vector(m_pIndex,u,Z,_out);
422}
423
424template <class T, enum MaterialOutput out>
425template <enum MaterialOutput _out>
426typename std::enable_if<_out==MaterialOutput::CauchyVectorN ||
427 _out==MaterialOutput::CauchyVectorM, gsMatrix<T>>::type
428gsMaterialMatrixIntegrateSingle<T,out>::eval3D_impl(const gsMatrix<T>& u, const gsMatrix<T>& Z) const
429{
430 return m_materialMat->eval3D_CauchyVector(m_pIndex,u,Z,_out);
431}
432
433
434template <class T, enum MaterialOutput out>
435template <enum MaterialOutput _out>
436typename std::enable_if<_out==MaterialOutput::PStress ||
437 _out==MaterialOutput::PStressN ||
438 _out==MaterialOutput::PStressM, gsMatrix<T>>::type
439gsMaterialMatrixIntegrateSingle<T,out>::eval3D_impl(const gsMatrix<T>& u, const gsMatrix<T>& Z) const
440{
441 return m_materialMat->eval3D_pstress(m_pIndex,u,Z,_out);
442}
443
444} // end namespace
gismo::gsFunctionSet
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition gsFunctionSet.h:219
gismo::gsGaussRule
Class that represents the (tensor) Gauss-Legendre quadrature rule.
Definition gsGaussRule.h:28
gismo::gsMaterialMatrixBase
This class defines the base class for material matrices.
Definition gsMaterialMatrixBase.h:33
gismo::gsMaterialMatrixIntegrateSingle
This class serves as the integrator of material matrices, based on gsMaterialMatrixBase.
Definition gsMaterialMatrixIntegrate.h:103
gismo::gsMaterialMatrixIntegrateSingle::eval3D_impl
std::enable_if< _out==MaterialOutput::VectorN||_out==MaterialOutput::VectorM, gsMatrix< T > >::type eval3D_impl(const gsMatrix< T > &u, const gsMatrix< T > &Z) const
Specialisation of eval3D for vectors.
Definition gsMaterialMatrixIntegrate.hpp:410
gismo::gsMaterialMatrixIntegrateSingle::eval_into_impl
std::enable_if< _out==MaterialOutput::Density, void >::type eval_into_impl(const gsMatrix< T > &u, gsMatrix< T > &result) const
Specialisation of eval_into for densities.
Definition gsMaterialMatrixIntegrate.hpp:170
gismo::gsMaterialMatrixIntegrateSingle::multiplyZ_into
void multiplyZ_into(const gsMatrix< T > &u, index_t moment, gsMatrix< T > &result) const
Multiplies the evaluation at z=0 with 2.0/(moment+1) * Thalf^(moment + 1)
Definition gsMaterialMatrixIntegrate.hpp:367
gismo::gsMaterialMatrixIntegrateSingle::eval_into
void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Implementation of eval_into, see gsFunction.
Definition gsMaterialMatrixIntegrate.hpp:160
gismo::gsMaterialMatrixIntegrateSingle::integrateZ_into
void integrateZ_into(const gsMatrix< T > &u, const index_t moment, gsMatrix< T > &result) const
Integrates through-thickness using Gauss integration.
Definition gsMaterialMatrixIntegrate.hpp:280
gismo::gsMaterialMatrixIntegrateSingle::multiplyLinZ_into
void multiplyLinZ_into(const gsMatrix< T > &u, const index_t moment, gsMatrix< T > &result) const
Uses the top and bottom parts of the thickness to compute the integral exactly.
Definition gsMaterialMatrixIntegrate.hpp:333
gismo::gsMaterialMatrixIntegrateSingle::_eval3D
gsMatrix< T > _eval3D(const gsMatrix< T > &u, const gsMatrix< T > &Z) const
Integrateuates the base class in 3D.
Definition gsMaterialMatrixIntegrate.hpp:401
gismo::gsMaterialMatrixIntegrateSingle::_eval
gsMatrix< T > _eval(const gsMatrix< T > &u) const
Integrateuates the base class in 3D, but on Z=0.
Definition gsMaterialMatrixIntegrate.hpp:393
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsQuadRule::mapTo
virtual void mapTo(const gsVector< T > &lower, const gsVector< T > &upper, gsMatrix< T > &nodes, gsVector< T > &weights) const
Maps quadrature rule (i.e., points and weights) from the reference domain to an element.
Definition gsQuadRule.h:177
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsGaussRule.h
Provides the Gauss-Legendre quadrature rule.
gsMaterialMatrixIntegrate.h
Provides an evaluator for material matrices for thin shells.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::freeAll
void freeAll(It begin, It end)
Frees all pointers in the range [begin end)
Definition gsMemory.h:312